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FOREWORD 


The  Transportation  and  Facilities  Research  Division  of  the  Agricultural 
Research  Service  conducts  research  to  find  ways  to  hold  down  the  costs 
of  physical  distribution  of  products  between  farms  and  consumers.  It 
seeks  to  determine  and  bring  about  the  adoption  of  the  most  efficient 
facilities,  equipment,  and  methods  for  moving  these  products  through 
distributive  channels. 

During  the  last  30  years,  as  grain  elevators  have  become  much  larger, 
i there  have  been  a number  of  structural  failures.  These  failures  appar- 
| ently  have  been  caused  by  excessive  pressures  on  the  walls  of  bins  or 
i silos  by  bulk-stored  grain.  Many  of  these  failures  have  resulted  in 
substantial  financial  losses.  Conversely,  there  is  reason  to  believe, 
because  of  inadequate  data  on  the  pressures  exerted  by  bulk-stored  grain 
under  dynamic  conditions,  that  many  elevators  constructed  in  recent 
years  have  a built-in  "safety  factor"  that  has  unnecessarily  increased 
construction  costs.  Most  elevator  design  engineers  still  compute  the 
structural  strengths  required  in  these  facilities  on  the  theories  of 
static  pressures  proposed  by  Janssen  in  1895  and  modified  by  Koenen  in 
1896.  In  1966,  research  was  initiated  under  a contract  with  Dr.  Joel  D. 
Isaacson,  St.  Louis,  Mo.,  to  devise  a method,  through  the  use  of  analytic 
and  algebraic  models,  for  computing  the  pressures  exerted  by  bulk-stored 
grain  under  dynamic  conditions.  This  publication  contains  the  contractor's 
report  as  it  was  written.  The  views  expressed  are  those  of  the  contractor 
and  do  not  necessarily  represent  the  views  of  the  Division  and  the 
Department  of  Agriculture. 


William  C.  Crow 
Director 

Transportation  and  Facilities 
Research  Division 
Agricultural  Research  Service 
U.S.  Department  of  Agriculture 
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INVESTIGATIONS  ON  THE  BASIC  THEORY  OF  STATIC  AND  DYNAMIC 
PRESSURE  PHENOMENA  IN  GRAIN  UNDER  CONDITIONS  OT  STORAGE 
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By 
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Joel  D.  Isaacson 


CHAPTER  1 


INTRODUCTION  AND  SUMMARY 


1.1.  INTRODUCTION 

The  determination  of  static  and  dynamic  pressures 
exerted  by  grain  on  the  containing  deep  bin  structure 
continued  to  be  a hazardous  task  for  the  designing  engineer. 
Mathematical  and  computer  investigations  on  the  basic 
theory  underlying  grain  pressure  mechanisms  are  presented 
in  this  report. 


1.1.1.  The  Nature  of  the  Problem 

Grain  pressures  in  deep  bins  were  first  calculated 
as  for  a semi-liquid  of  the  same  density  as  the  grain.  Thus 
for  example,  the  lateral  pressure  was  calculated  from  the 
hydrostatic  formula. 


p(y>  = yy 


(l.D 


where, 

p (y)  - the  lateral  pressure  exerted  by  the  grain 
on  the  wall  at  the  depth  y 
Y - the  density  of  the  grain 

- the  depth  of  the  cross-section  calculated. 


y 


This  method  was  inadequate  because  many  structures 
buckled  under  the  vertical  load  arising  from  the  friction 
of  the  grain  on  the  walls  - a phenomenon  unknown  at  the 
time.  On  the  other  hand,  it  was  realized  later  that  the 
calculations  based  on  hydrostatic  pressure  distribution 
yielded  exaggerated  results  compared  with  the  actual 
lateral  pressure* 

The  developments  in  the  field  of  soil  mechanics 
for  calculating  earth  pressure  gave  rise  to  the  definition 
of  a new  factor,  kQ,  "the  coefficient  of  earth  pressure 
at  rest",  which  was  borrowed  for  use  in  the  grain  pressure 
problem.  Equation  (1. 1)  thup  became. 


p (y ) = kQYy 


(1.2) 


where  kQ  was  usually  determined  according  to  Rankine  as 


l-sin8 

l+sin0 


(1.3) 


where  0 is  the  angle  of  internal  friction  of  the  grain. 

This  improved  the  estimation  of  the  lateral  pressure,  but 
still  did  not  take  into  account  the  vertical  load  on  walls. 
The  whole  weight  of  the  grain  was  erroneously  assumed  to 
be  transmitted  to  the  bottom. 

The  next  development  was  proposed  by  Janssen  (1895) 
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took  into  account  the  grain-wall  friction. 


P(y)  = [l-exp{- (-g^)y}]  (1.4) 

where , 

R - the  "hydraulic  radius"  of  the  bin  cross-section 

y - the  coefficient  of  friction  of  the  grain  on 
the  wall 

k - Rankine  coefficient  of  earth  pressure  at 'rest, 
proposed  for  this  case  by  Koenen  (1896) 

Since  that  time  a considerable  number  of  formulas 
have  been  proposed  by  various  investigators.  Most  of 
them  give  similar  results  to  that  of  the  Janssen-Koenen 
method,  and  consequently  the  latter  has  remained  in  widest 
use. 

Over  a long  period  it  appeared  that  the  problem 
has  been  solved;  however,  during  the  last  three  decades, 
when  grain  elevators  became  much  larger  and  more  numerous, 
some  failed  due  to  excessive  pressures.  Investigations 
have  shown  that  during  discharge  operations  so  called  "dynamic 
pressure"  develop  that  may  be  twice  as  high  as  the  calculated 
Janssen's  pressure.  See  Figure  1.1. 

. Whereas  the  conventional  formulas  are  more  or  less 
applicable  to  static  conditions,  there  do  not  exist 
satisfactory  methods  of  estimating  grain  pressures  under 
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1 - Hydrostatic  pressure  curve,  Eq.  (1*1) 

2 - Earth  pressure  curve,  Eq.  (1.2) 

3 - Conventional  pressure  curve  (Janssen's),  Eq.  (1*4) 

4 - Typical  empirical  dynamic  pressure  curve 


Lateral  Pressure  p 


Figure  1. 1 Geometrical  representation  of  the  nature  of 
the  problem. 
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Jlynamic  conditions.  An  attempt  to  rectify  this  situation 
is  made  in  this  report. 


1.2  SUMMARY 

The  following  is  a brief  summary  of  topics  treated 
in  this  report. 

Chapter  2 is  introductory,  giving  general  definitions 
of  factors  associated  with  grain  pressure  systems.  It  also 
defined  deep  bin  in  quantitative  terms,  and  devises  a 
compact  classification  of  granular  materials. 

In  chapter  3,  the  basic  one-dimensional  pressure 
mechanism  is  formulated,  first  in  analytic  and  then  in 
algebraic  terms.  The  analytic  model  is  based  on  the 
solution  of  an  ordinary  linear  differential  equation  of  the 
first  order  with  (possibly)  variable  coefficients.  A 
practical  solution  of  the  analytic  model  is  given,  accompanied 
by  a computer  program. 

The  algebraic  model  is  given  in  a compact  form, 
involving  series  of  determinants  and  matrices,  readily 
applicable  to  numerical  calculations.  The  models  are  very 
general,  allowing  both  static  and  dynamic  conditions,  varied 
bin  geometries,  and  surcharge/no  surcharge  conditions.  They 
are  compared  for  various  cases.  For  the  case  of  constant 
coefficients  the  analytic  and  the  algebraic  models  are 
shown  to  be  mutually  identical  as  well  as  identical  with 
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Janssen* s solution.  For  variable  coefficients  (a  linearly 
increasing  density-function,  a hyperbolically  decreasing 
ratio-function,  and  a constant  wall-friction  parameter) 
the  models  stay  identical.  A special  choice  of  the 
parameters  leads  to  complete  identity  with  the  Reimberts ' 
semi-empirical  solution.  The  models  are  shown  to  be  proto- 
types of  many  solutions,  from  multi-dimensional  dynamic  to 
the  one-dimensional  static,  the  Reimberts',  Janssen's, 
Rankine's  and  the  hydrostatic  solutions. 

The  characteristic  functions  are  discussed  in 
chapter  4.  These  are  the  density,  ratio  and  friction 
functions,  giving  the  unit  weight,  the  ratio  between  the 
lateral  to  the  vertical  pressures,  and  the  ratio  between 
the  frictional  stress  and  the  lateral  pressure,  as  functions 
of  the  depth  coordinate.  Typical  characteristic  functions 
associated  with  static  and  dynamic  conditions  are  derived 
and  analyzed.  Some  practical  recommendations  follow. 

Numerical  and  computer  models  are  given  in  chapter  5. 
Three  examples  of  numerical  models  based  on  the  analytic 
model  are  given  first.  The  models  are  complete  with  computer 
programs  and  printed  results.  Next  are  given  four  examples 
of  computerized  graphic  representations  of  various  analytic 
models.  Lastly,  six  examples  of  digital  simulations  based 
on  the  algebraic  model.  These  models  are  also  complete 
with  computer  programs  and  printed  results. 
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Chapter  6 outlines  methods  for  utilizing  the  basic 


one-dimensional  models  to  solve  models  of  higher  complexities. 
The  first  section  gives  formal  generalizations  of  the  one- 
dimensional models  and  characteristic  functions.  In 
section  6.2  a computer-oriented  method  of  solving  two- 
dimensional  problems  by  iterative  use  of  one-dimensional 
models  is  introduced.  Section  6.3  discusses  the  complexities 
introduced  into  the  grain  pressure  system  due  to  feedback 
and  coupling  mechanisms.  The  problem  is  handled  through 
graph- theoretic  techniques.  Problems  of  this  category  can 
be  represented  by  directed  graphs  and  solved  by  controlled 
iteration  of  one-dimensional  problems. 

Chapter  7 introduces  an  approach  based  on  topological 
considerations.  A brief  presentation  of  the  principles 
leads  to  intermediate  results  which  show  that  the  maximum 
dynamic  lateral  pressure  in  typical  circular  cylindrical 
grain  bins  should  occur  at  about  one-third  of  the  total 
height  H above  the  bottom  and  decrease  sharply  toward  a 
point  about  0.12  H above  the  bottom.  The  last  section 
describes  briefly  a computer  system  developed  under  this 
contract  to  produce  grain  pressure  simulations  in  the  form 
of  16mm  motion-pictures. 

Appendix  A gives  an  exhaustive  bibliography  on 
grain  pressure.  The  literature  survey  covers  material 
published  since  1883  on  a world-wide  basis.  Among  the 
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countries  represented  in  the  survey  are  (in  alphabetical 
order):  Algeria,  Argentina,  Austria,  Canada,  England, 

France,  Germany,  Hungary,  Ireland,  Israel,  Italy,  Mexico, 
the  Netherlands,  Poland,  South  Africa,  the  Soviet  Union, 
Sweden  and  the  U.S. 

Appendix  B gives  a complete  program  listing  of 
a package  for  automatic  three-dimensional  plots. 
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CHAPTER  2 


DEFINITIONS 

2.1  GENERAL  DEFINITIONS 

(1)  A bin  is  a container  that  bounds  geometric  figures  belonging 
to  the  set  of  all  right  cylinders  that  satisfy  either  of  the  following 
conditions: 

a.  The  directrix  of  the  cylindrical  surface  is  an  arbitrary 
closed  convex  curve. 

b.  The  directrix  is  a pair  of  infinite  parallel  straight 
lines;  the  bin  is  then  called  an  infinite  bin. 

The  cylindrical  surface  is  called  the  wall( s ) or  boundary  of  the  bin. 

The  lower  base  is  called  the  (flat)  bottom. 

The  upper  base  is  called  the  top,  and  physically  may  or  may  not  exist. 

(2)  The  hydraulic  radius  of  a bin  is  in  the  respective  cases 

above: 

a.  The  ratio  of  the  area  bounded  by  the  closed  directrix 
to  the  length  of  the  directrix. 

b.  The  ratio  of  the  area  bounded  by  the  rectangular  cross- 
section  of  a finite  segment  of  an  infinite  bin,  to  the  total  length  of 
the  walls  included  in  the  segment. 

(3)  A bin  is  shallow  or  deep  according  to  whether  the  ratio  of 
its  height  to  the  hydraulic  radius  is  less  than  or  equal  to  or  larger  than 
some  arbitrary  positive  real  number,  respectively. 

(U)  A bin  is  called  a grain  bin  if  its  interior  is  partially  or 
completely  occupied  by  grain.  (See  def.  6.) 


9 


(5)  Granular  material  is  a mass  consisting  of  a collection  of 


solid  particles  and  has  the  following  properties: 

a.  A fixed  natural  angle  of  repose.  (See  def.  7.) 

b.  If  bounded  by  a deep  bin,  it  maintains  gravity  flow 
when  some  region  of  the  bottom  is  removed. 

c.  If  bounded  by  a deep  bin,  it  transforms  to  a geometric 
figure  of  constant  slopes  when  the  walls  are  gradually  removed. 

(6)  Grain  is  a collection  of  plant  seeds  that  satisfies  the 
properties  of  granular  materials.  A s-ngle  member  of  the  collection  is 
called  a kernel . 

(7)  The  natural  angle  of  repose  of  a granular  material  is  the 
angle  of  the  steepest  slope  of  the  right  circular  cone  with  horizontal 
base  that  the  granular  material  may  form  under  gravity. 

(8)  The  lateral  pressure  at  a point  on  the  wall  surface  is  the 
sum  of  the  horizontal  forces  exerted  by  the  grain  on  a square  neighbor- 
hood of  this  point , of  unit  area. 

(9)  The  frictional  stress  at  a point  on  the  wall  surface  is  the 
sum  of  the  vertical  forces  exerted  by  the  grain  (due  to  lateral  pressure 
at  this  point)  on  a square  neighborhood  of  this  point,  of  unit  area. 

(10)  The  wall  friction  parameter  at  a point  on  the  wall  surface 
at  a given  instant,  is  the  ratio  between  the  frictional  stress  and  the 
lateral  pressure  at  this  point,  at  the  given  instant. 

(11)  The  vertical  pressure  at  a point  on  a horizontal  cross-section 
of  the  grain  is  the  sum  of  the  vertical  downward  forces  exerted  by  the 
grain  on  a square  neighborhood  of  this  point,  of  unit  area. 

(12)  The  bottom  pressure  is  the  vertical  pressure  on  the  bottom. 
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2 


DEFINITION  OF  A DEEP  BIN 


The  wall  friction  effect  under  static  conditions  characterizes  the 
main  difference  "between  deep  and  shallow  "bins.  It  is  therefore  desirable 
to  use  the  wall  friction  effect  as  a measure  of  the  "deepness"  of  a "bin. 

For  this  purpose  Janssen* s equation  may  "be  used. 

The  asymptotic  behavior  of  the  pressure  is  due  to  the  exponential 
function 

Pk 

£(y)  = exp{  -( — 2)  y}  (2.1) 

R 

where  0 < ( y ) <_  1 for  0<  y < 00 . 


By  definition,  a bin  is  deep  if  its  total  depth  H satisfies  the 

inequality  H > H , where  H is  a depth  such  that  4(H)  = 0.05* 

d d d 

In  other  words,  a bin  can  be  considered  to  be  deep  if  its  depth  is 

sufficiently  large  such  that  at  least  95  percent  of  the  wall  friction 

effect  is  developed  at  some  level  above  the  bottom,  and  any  load  below  this 

level  is  essentially  carried  by  the  walls  alone,  due  to  friction. 

To  solve  for  the  ratio  H /R,  insert  in  Eq.(2.l), 

d u 

pk 

C(H  ) = exp{  -{—)  H } = 0.05 
d R <3 

and  thus 


(2.2) 


For  example,  in  terms  of  the  ratio  H /D  for  circular  cylindrical 

d 


bins,  one  obtains 


H 


= - x — 

yk 


(2.3) 


11 


Typical  values  of  k and  y for  grain  are  respectively  0:333  and  0.U50, 

o 

and  therefore 

Hd  3 1 

= - x = 5, 

D It  0.U50x0.333 

so  that  for  practical  considerations  a circular  cylindrical  grain  bin 
is  deep  if 

H 

- 1 5.  ( 2.k ) 

' D 

2.3  CLASSIFICATION  OF  MATERIALS 

Granular  materials  are  sometimes  referred  to  as  semi-fluids  or 
bulk-solids , probably  to  suggest  their  intermediate  properties.  It  seems 
that  granular  materials  subject  to  gravity  pressures  indeed  possess 
properties  whose  ranges  lie  between  those  of  fluids  and  monolithic  solids. 
Table  2.1  presents  a classification  of  materials  with  respect  to  gravity 
pressures  in  storage  containers.  From  this  classification  there  arises 
a compact  formulative  way  to  designate  materials  by  their  gravity-pressure 
properties.  For  example,  a material  satisfying  the  hydrostatic  equation 
for  lateral  pressure  is  designated  by  I-A1-B1-C1-D1-E1 , i.e.  water.  This 
kind  of  designation  is  called  a material-type . More  than  one  material  is 
associated  with  a material -type.  In  general,  higher  order  characters  in  a 
material-type  indicate  more  complex  materials. 

Table  2.2  presents  some  examples  of  material-types. 
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TABLE  2.1  Classification  of  materials  with  respect  to  gravity  pressures  in  storage  containers 


RANGE  OF  VALUES 

W 

(M 

var. 

X X 

nJ  nj 

6 g 

V'  Vi 

' 7*- 

• V*  v.  ; 

i o o' 

const. 

O O O O 

r r fr  r 

Q 

O 

var. 

o<c 

H 

const. 

c — o 
c=o 

o<c 

O 

CD 

c* 

var. 

t=|r\i  t=|cvj 

i V Vi 

to  CD  1 

1 V Vi 

0 O 

H 

const. 

0 V V 

ll  © © Vi 

CD  V V t=1fM 

0 O 

PQ 

CvJ 

var. 

pH  pH 

i Vi  Vi  i 

t ^ r<  1 

i Vi  Vi  i 

O O 

pH 

const. 

pH  pH 

0 Vi  Vi  0 

ll  X ^ II 

^ Vi  Vi  ^ 

O O 

< 

(\j 

var. 

o<k<l 

°<k<l 

r-H 

const. 

k=l 

o<k<l 

°<k<l 

k=o 

MATERIAL 

I fluids 

II  non- 
cohesive 
granular 
materials 

III  cohesive 
granular 
materials 

IV  monolithic 
solids 

4) 

(0 

O 

cu 

4) 

U 

o 

D 

00 

g 

nJ 


cJ 

u 

D 

4-> 

C 

i 

CD 


4) 

M 

O 

•5 


I I 
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ratio  between  lateral  and  vertical  pressures  at  any  level  c - cohesion 

ratio  between  frictional  stress  and  lateral  pressure  at  any  level  y - unit  weight 


TABLE  2.2  Examples  of  material-types  of  materials  stored  in  bins 


Material-type 

Material 

a 

I-A1-B1-C1-D1-E1 

water,  oil 

b 

I I-A1-B1-C1-D1-E1 

dry  sand,  aggregate,  Janssen* s grain 

c 

II-A2-B2-C2-D1-E2 

grain 

d 

III-A2-B2-C2-D1-E1 

Portland  cement,  flour 

e 

III-A2-B2-C2-D2-E2 

silage 

f 

IV-A1-B1-C1-D0-E1 

conglomerated  grain,  hardened  concrete  1/ 

1/  Materials  may  change  material- types  due  to  chemical,  organic, 
biological  and  environmental  factors.  An  extreme  example  is  concrete. 

When  poured  into  the  forms  it  is  a mixture  of  components  (water,  sand, 
aggregate  and  cement)  belonging  to  material-types  (a),  (b),  (b),  (d), 
respectively.  Nevertheless,  the  first  component  is  qualitatively  dominant 
and  poured  concrete  is  closest  to  material -type  (a).  However,  after  a 
relatively  short  period  the  concrete  hardens  and  gradually  becomes  material- 
type  (f). 

Similar  changes,  although  less  critical,  may  occur  in  grain,  flour 
or  silage,  under  certain  conditions. 

The  concerns  of  this  study  are  grain-types:  II-A2-B2-C2-D1-E2  and 
simpler  types. 
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CHAPTER  3 


BASIC  PRESSURE  MECHANISM 

3.1  ANALYTIC  MODEL 

3.1.1  Introduction 

For  clarity,  the  one-dimensional  problem  is  dealt 
with  first.  In  the  one-dimensional  problem  all  variables 
are  considered  functions  of  the  depth  coordinate  - y - alone. 
For  simplicity,  the  infinite  bin  case  is  considered. 

Similar  treatment  can  be  applied  to  any  other  bin  with  a 
symmetrical  cross-section,  such  as  a circle,  a square,  or 
a regular  polygon. 

3.1.2  Derivation  of  the  Governing  Differential  Equation 
Consider  a horizontal  differential  slice  of  height 

dy  at  depth  y in  a unit  segment  of  an  infinite  bin,  having 
width  2p.  The  vertical  pressure  acting  on  the  top  surface 
of  the  prism  is  q,  the  lateral  pressure  exerted  by  the 
prism  of  grain  on  the  walls  is  p,  and  the  resulting 
frictional  stress  is  s.  Let  dq  denote  the  increase  in 
vertical  pressure  along  the  interval  (y,  y + dy) . Then 
the  increase  in  the  resultant  of  the  vertical  forces  acting 
on  the  prism  in  a downward  direction,  along  this  interval, 
is 

A dq  = y A dy  - s L dy  (3.1) 
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where 


A - cross-sectional  area 

L - that  portion  of  the  cross-sectional  perimeter 
which  is  bounded  by  the  walls 
y - unit  weight  of  grain 
Rearrangement  of  Eg.  (3.1)  yields 


Define 


and  thereby 


Insert  Eg. 


dq  L 

d?  + S A = Y- 

(3.2) 

k = H 

(3.3) 

q 

* = 1 

(3.4) 

p 

CM 

II 

P4 

(3.5) 

s = X k q. 

(3.6) 

(3.6)  into  Eg.  (3.2)  to  obtain 


g + 1 X k q = Y (3.7) 

and  assume  that  R,  k,  X,  y are  constants  or  functions  of  y 
alone . Denote 

e = 0.8) 

Then,  the  most  general  eguation  is  of  the  form 


+ 3 (y ) g = y (y). 


(3.9) 
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— H^q.  (3.9)  is  an  ordinary  linear  differential  equation  of 
the  first  order,  where  y and  q are  the  independent  and 
dependent  variables  respectively,  and  3 and  y are  functions 
of  y alone,  or  constants. 


3.1.3  Initial  Conditions  and  Solutions 


To  find  the  general  solution  -of  Eq.  (3.9)  consider 
first  the  homogeneous  linear  equation 

g + eq  = o. 

Its  variables  are  separable,  thus 

^2.  + gdy  = 0 

q * 


and  the  solution  is 


1 


q = c exp (-/  3 dy) 

where  c is  a constant. 

Now  substitute  in  the  non-homogeneous  equation,  the 
expression 

q = w exp {-J$dy) 

in  which  w,  a function  y,  replaced  the  constant  c.  The 
equation  becomes 


g Bxpi-Jt 


$dy)  = y. 
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whence 


w = C + J y exp  (fgdy)  dy. 
The  general  solution  is  therefore 
q=C  exp(-f3dy)  + exp  (- 


-f BdyjyC  exp 


3dy)  dy  (3.10) 


where  C is  an  arbitrary  constant.  In  order  to  determine 
C uniquely,  an  initial  condition  must  be  specified.  In 
this  case  (a  one-dimensional  problem)-,  it  is  necessary 
and  sufficient  to  specify  at  the  initial  point  y = 0, 
the  following  initial  conditions: 

0 without  surcharge 


g (o ) = 


( 

( l/2yol 


(3.11) 


with  triangular 

surcharge  of  height  h? 

Y is  constant, 
o 


(1)  Solution  without  surcharge.  Imposing  the 
first  initial  condition  on  Eq.  (3.10)  one  obtains 


C = 


-J y exp  ( J t 


Y exp  ( j 3 dy)  dy 


(3.12) 


y=o 


and  by  Eqs.  (3.3),  (3.10),  and  (3.12),  a complete  solution 

for  the  lateral  pressure  is  of  the  form 


p=k  exp (-/ 3 dy) 


- /*3  dy)|^Yexp  (J&  dy)  dy  -J y exp  ( Jb  dy ) dyj^ 


(3.13) 


(2)  Solution  with  surcharge.  Imposing  the  second 
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initial  condition  on  Eq.  (3.10)  one  obtains 


C = { l/2yQh  exp ( |B  dy)  -jy  exp  ( 13  dy)  dy  }J 


(3.14) 


y=o 


and  a complete  solution  is  of  the  form 
p = k exp (-J 3 dy) 


y exp ( 3j dy)  dy  + {l/2YQh  exp ( / 3 dy)- 


'/6 


-/ 


- Y exp  ( / 6 dy) 


dy}  "1  . 

y=oJ, 


(3.15) 


Notes . By  the  broad  definition  of  L and  A the 
solution  is  in  fact  applicable  not  only  to  the  infinite 
bin  case,  but  also  to  any  symmetrical  bin.  Since  case  (1) 
is  less  complicated,  unless  otherwise  indicated  Eq.  (3.13) 
will  be  hereinafter  referred  to  as  the  analytic  model. 

Eq.  (3.13)  gives  a complete  analytic  solution  to 
the  one-dimensional  problem,  provided  3 and  y are  known 
and  can  be  expressed  in  rather  simple  forms.  It  will  be 
shown  later  (see  Chapter  4*  ) that  plausible  assumptions 
can  be  made  regarding  the  forms  of  these  functions.  How- 
ever, in  general,  it  is  best  to  solve  Eq.  (3.13)  numerically 
using  a procedure  that  admits  arbitrary  forms  of  3 and  y. 

A procedure  such  as  this  is  presented  in  the  next  section. 
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3.2  NUMERICAL  SOLUTION  (RUNGE-KUTTA)  FOR  THE  ANALYTIC 

MODEL 

3.2.1.  Introduction 

Initial  value  problems  can  be  stated  as  follows: 
Given 


l 


I.C.:  y(xo)=yo. 


(3.16) 


(3.17) 


find  a solution  y = g (x) , such  that  yQ=g (xq) . (Here, 
x and  y are  the  independent  and  dependent  variables, 
respectively) . 

To  express  the  grain  pressure  problem  in  these 
terms  it  is  necessary  to  rewrite  Eqs.  (3.9)  and  (3.11) 
as  follows. 


(^3-  = y (y)  - B (y)  q (3. 18) 

q ( 0 ) = constant  (3.19) 

The  R.H.S.  of  Eq.  (3.18)  can  be  any  arbitrary  function 
of  y and  q,  thus  allowing  a great  flexibility  in  choosing 
3 and  y.  The  constant  in  the  initial  condition  depends 
on  the  surcharge.  It  is  usually  zero. 
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With  the  advent  of  modern  high-speed  digital 
computers,  numerical  methods  for  the  solution  of  initial 
value  problems  in  ordinary  differential  equations  have 
found  widespread  application.  A variety  of  methods  and 
special-purpose  simulation  languages  are  available. 

Benyon's  article  (1)  on  numerical  methods  for  digital 
simulation  reviews  many  of  these  methods.  1/  Further 
information  can  be  found  in  references  (2^)  through  (4_2)  . 

After  considerable  investigation  it  was  found  that  for 
the  purpose  of  solving  the  analytic  model,  a Runge-Kutta 
method  is  well-suited.  A brief  description  of  the 
method  follows. 

3.2.2  A Fourth-Order  Runge-Kutta  Method 

Runge-Kutta  methods  involve  essentially  replacing 
a truncated  Taylor* s series  expansion  of  the  solution, 
given  in  terms  of  derivatives,  by  an  approximation  in 
terms  of  f(x,y)  only.  The  derivation  of  the  method  can 
be  found  in  many  books  on  numerical  analysis,  for  example  ( 24 ). 

Using  a fourth-order  Runge-Kutta  integration 
procedure  the  ordinary  differential  equation,  dy/dx=f(x,y) 
with  initial  condition  y(xQ)  = yQ  is  solved  numerically. 

This  is  a single-step  method  in  which  the  value  of 
y at  x = xn  is  used  to  compute  = y(xn+-^)  anc^  earlier 

values  y , , y etc.  are  not  used. 

Jn-1'  Jn-2 

17 Underscored  numbers  in  parentheses  refer  to 
References,  p.  150. 
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The  relevant  formulas  are: 

yn+1  = yn  + 1/6  (kQ  + 2kx  + 2k2  + k3)  (3.20) 

Where , for  step  size  h: 

k = hf (x  , y ) 
o nfjrn 

kx  = hf(xn  + h/2,  yn  + kQ/2) 

(3.21) 

k2  = hf(xn  + h/2,  yn  + k±/2) 

k3  = hf(xn  + h,  yn  + k2) 

There  formulas  were  programmed  and  implemented 
on  the  computer.  A subroutine  RUNKUT  follows. 
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oooooooooooooooooooooooooo 


SUBROUTINE  RUNKUT 


PURPOSE 

INTEGRATES  THE  FIRST  ORDER  DIFFERENTIAL  EQUATION  OF  THE 
MODEL,  DY/DX=FUN(X,Y)  AND  PRODUCES  A TABLE  OF  INTEGRATED  VALUES 


USAGE 

CALL  RUNKUT1FUN,H,XI ,YI,K,N»VEC) 

DESCRIPTION  OF  PARAMETERS 

FUN-SUPPLIED  FUNCTION  SUBPROGRAM  WITH  ARGUMENTS  X,Y 
(DEPTH  AND  LATERAL  PRESSURE)  WHICH  GIVES  DY/DX 
H -STEP  SIZE  ALONG  DEPTH  COORDINATE 
XI  -INITIAL  VALUE  OF  X (ZERO  USUALLY) 

Y I -INITIAL  VALUE  OF  Y WHERE  YI=Y(XI)  (USUALLY  ZERO, 

OR  DEPENDS  ON  THE  SURCHARGE) 

K -THE  INTERVAL  AT  WHICH  THE  COMPUTED  VALUES  ARE  TO  BE  STORED 
N -THE  NUMBER  OF  VALUES  TO  BE  STORED 

VEC-THE  RESULTANT  VECTOR  OF  LENGTH  N IN  WHICH  COMPUTED 
VALUES  OF  Y ARE  TO  BE  STORED 


SUBROUTINE  RUNKUT ( FUN , H , X I ,YI ,K,N,VEC) 
DIMENSION  VEC(I) 

H2=H/2. 

Y= Y I 
X=X  I 

DO  2 1=1, N 
DO  1 J=1 , K 
T1  = H*FUN<  X,Y) 

T2=H*FUN( X+H2,Y  + Tl/2*  ) 

T3=H*FUN(  X+H2,  Y+T2/2.  ) 

T4=H*FUN(  X+H , Y + T3 ) 

1 X=X+H 

2 VEC ( 1 ) =Y 

Y=  Y+(T1+2.*T2+2.*T3+T4) /6. 

RETURN 

END 
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3.3 


ALGEBRAIC  MODEL 


In  addition  to  the  analytic  model,  it  seems 
desirable  to  supply  an  alternative  algebraic  solution 
that  is  free  from  any  analytical  difficulties  that 
may  normally  be  associated  with  a complicated  expression 
such  as  Eq.  (3.13).  Another  reason  for  seeking  an 
algebraic  solution  is  to  make  practical  calculations 
more  readily  applicable  to  programming  and  processing 
by  digital  computers. 

3.2.1  Construction  of  the  Basic  Model 

Consider,  again,  a unit  segment  of  the  infinite 
bin  having  width  and  height  2p  and  H,  respectively.  For 
reasons  of  symmetry,  treat  one-half  only.  Divide  H into 
u slices  each  of  height  Ay.  Starting  from  the  top, 
associate  with  each  slice  an  ordinal  number  i,  (l<i<u,i 
is  an  integer),  and  local  values  for  various  parameters, 
such  as:  unit  weight  , wall-friction  parameter  A^, 
and  ratio  parameter  K^. 

NOMENCLATURE: 

y^  - average  unit  weight  of  the  i-th  slice 
A - cross-sectional  area,  A=l.  p=  p 

- weight  of  the  i-th  slice,  G^  = y^AAy 
Q.  - resultant  of  the  downward  vertical  forces 

l 

acting  on  the  i-th  slice 
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- resultant  of  the  lateral  forces  acting 
on  the  i-th  slice 

- resultant  of  the  vertical  frictional 

forces  exerted  by  the  i-th  slice  on  the 

wall  due  to  F. 

1 

K.  - ratio  between  F.  and  Q.,  K.  = F./Q. 

1 1 l i i 

A^  - wall  friction  parameter  at  the  i-th  slice 

(ratio  between  S.  and  F.,  A.  = S./F.) 

1 l'  1 a/  1 

- average  vertical  pressure  on  a horizontal 
cross-section  passing  through  the  center 
of  gravity  of  the  i-th  slice 

p^  - average  lateral  pressure  exerted  by  the 
i-th  slice  on  the  wall 
s ^ - average  frictional  stress  at  the  wall, 
exerted  by  the  i-th  slice 
k.  - ratio  between  p.  and  q.,  k.  = p./cr- 

- accumulated  load  per  unit  length  carried 
by  the  wall  at  the  i-th  slice 

p - hydraulic  radius 

Isolate  slice  #1  and  examine  the  forces  acting 

upon  it: 

(i)  weight  of  slice  - (downwards) 

(ii)  resultant  of  lateral  forces  - fi.=K1G1  (inwards) 

(iii)  resultant  of  frictional  forces  - 

(upwards  at  the  wall) 
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The  residual  resultant  of  downward  forces  is  therefore 

R1  = 

is  that  portion  of  that  is  transmitted  to  slice  #2, 
the  remaining  portion  being  carried  by  the  wall.  The 


same 

process 

can 

be 

applied 

to  each  slice 

i as  follows: 

slice 

# 

G. 

l 

K. 

l 

X. 

i 

F'. 

l 

S! 

i 

R. 

i 

1 

Gi 

Ki 

KiGi 

xiKiGi 

G1(l-X1K1) 

2 

G2 

K2 

X2 

K2G2 

^2K2G2 

G2(1-X2K2) 

9 

• 

• 

9 

• 

• 

• 

9 

• 

• 

9 

• 

• 

• 

• 

• 

• 

9 

• 

• 

• 

i 

Gi 

K. 

i 

X. 

1 

K.G. 

l i 

X.K.G. 

ill 

G. (1-X.K. ) 

l ll 

• 

• 

• 

• 

• 

• 

• 

• 

9 

• 

• 

• 

• 

• 

u 

Gu 

K 

u 

X 

u 

K G 
u u 

X K G 
u u u 

G (1-X  K ) 
u u u 

where 

F! 

l 

and 

S! 

l 

are 

those  portions  of  F^ 

and  S^,  respectively 

due  to  alone. 

A portion  of  the  weight  of  each  slice  is  transmitted 
to  the  wall  by  friction.  The  residual  portion  is  trans- 
mitted downwards  to  the  successive  slice.  After  is 
transmitted  to  the  i+1  slice,  it  is  again  divided  into 
two  portions:  one  is  transmitted  to  the  wall  and  the  residual 
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is  transmitted  to  the  i+2  slice,  and  so  on.  Thus,  the 
complete  process  is  illustrated  as  follows:  2/ 


Slice  # 

1.  Slice  #1  - 

2.  Slice  #2  - 
from  slice  #1  - 

3.  Slice  #3  - 
from  slice  #2 


Vertical  Downward  Forces 


(G1(l-X1K1) 

f3 

°3  “jG2(1-X2K2) 

Igi(1-X1K1)  (1-X2K2) 


i.  Slice  #i  - 

from  slice  #i-l  - 
from  slice  #i-2  - 


from  slice  #1  - 


1Gi-l(1_Xi-lKi-l) 

'i-2  (1"Ai-2Ki-2)  (1_Xi-lKi-l) 


Qi  ^ 


<G1(1-X1K1)  (1-X2K2)  . . . (l-X^K.^) 


2/  The  derivation  is  made  for  a case  without  surcharge. 
Surcharge  may  be  simply  added  to  G^. 
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The  sum  of  the  vertical  forces  acting  (downwards) 


on  the  n-th  slice  is  therefore 

Q =Ayp  [y  +y  -,(1-1  J ) +y  ~(1-X  ~K  0)(1-X  ,K  ,)+... 
n -*Kt,n  'n-1  n- 1 n-1  'n-2  n-2  n-2  n-1  n-1 

. . .+y1(l-X1K1) (1-X2K2) . . . (l-Xn-1Kn_1) ] . (3.22) 

The  series  in  Eq.  (3.22)  characterizes  the  model 
sought.  It  is  desirable  to  find  a more  compact  form  for 
this  series  so  that  it  can  be  handled  and  analyzed 
conveniently.  The  series  can  be  represented  as  a series 
of  determinants  as  follows: 


y I 1 1 +y  -i 
’n1  1 'n-1 


0 


0 (1-X  , K ,) 

n-1  n-1 


+Y 


n-2 


0 (1-X  , K ,) 

n-1  n-1 


0 


0 (1-X  0K  0) 
n-2  n-2 


1 0 

1-X  n K V J 

n-1  n-1 


1-X  0K  0 
n-2  n-2 


. . . + y 


0 
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The  determinants  above  can  be  recognized  as  being 
associated  with  the  sequence  of  n submatrices  produced  from 
the  n-by-n  diagonal  matrix 


i 

1 1 0 . 0 
I 1 


by  successively  taking  the  1-by-l  upper  left  submatrix, 
2-by-2  upper  left  submatrix,  and  so  on.  Let  the  expression 
"det  Sub.  A"  denote  "the  determinant  of  the  i-th  submatrix 

l 

of  an  n-by-n  arbitrary  diagonal  matrix  A,  produced  by 
deleting  from  A all  rows  and  columns  whose  ordinal  number 
is  larger  than  i."  Then  the  series  of  determinants  can 
be  written  in  the  form  of  a compact  sum  as  follows: 
n 

I y ,,  -det  Sub.  [I-AK]  (3.23) 

' n+l-i  l 

where  I is  an  n-by-n  identity  matrix,  and  AK  is  a diagonal 
matrix  of  the  form 
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(3.24) 


AK=Diagonal  [0 f Xn-lKn-l' Xn-2Kn-2 ' * # * • 

Thus  Eq.  (3.22)  becomes 
n 

Qn=Ayp  E Yn+1_i  det  Subi  [I  - AK]  (3.25) 

1=1 


and  other  quantities  are  readily  obtainable  in  terms  of 
Eq.  (3.25)  as  follows: 

(i)  Average  vertical  pressure  at  the  n-th  slice 


= Ay  E Yn+2.-i^et  Sub  = [I  - AK]  . (3.26) 


(ii)  Average  lateral  pressure  exerted  on  the  wall  at 
the  n-th  slice 


p = k q ; 
^n  n^n 


, K p 
k _ nK 

n Ay 


(3.27) 


n 

= Knp  E ^n+l-idet  Subi[I  " XK] * 
i=l 


(3.28) 


(iii)  Average  frictional  stress  at  the  wall  at  the 
n-th  slice 


s = A p 
n n^n 


n 


= A K p E y , , .det  Sub.  I - AK] . 3.29 

n nK  . , ' n+l-i  l 

i=l 
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(iv)  Accumulated  load  per  unit  length  carried  by 


the  wall  at  the  n-th  slice 


n j 

= Z X .K.  Z y . , , .det  Sub. [I  - XK] . (3.30) 

j = l =>  3 i=l  =>  + 1-1  1 
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COMPARISON  BETWEEN  THE  MODELS  AND  RELATIONSHIPS 


3.  4 


TO  OTHERS 


Clearly  the  two  models  were  mathematically  derived 
differently,  though  both  represent  the  same  physical 
situation.  It  is  of  interest  to  carry  out  a preliminary 
investigation,  comparing  between  the  models,  as  well  as 
with  known  models  such  as  Janssen's  and  the  Reinberts*. 


3.4.I  Constant  Coefficients 

Start  out  with  the  simplest  case  where  all 
parameters  are  constant.  Let  X = y,  k,  y,  and  p be  constants 
Also,  let  all  the  non-zero  elements  of  the  diagonal  matrix 
XK  be  the  same,  such  that  (XK)^=yK,  where  K=—  Ay  is  constant 
All  the  derivations  are  with  respect  to  lateral 
pressure  only. 

(1)  Analytic  solution.  The  lateral  pressure  p 
is  given  by  Eq.  (3.13).  Substitute  the  constant  coefficients 
into  Eq.  (3.13)  to  obtain 


where  b 


is 


a constant. 


Then 


p = ky  exp (-by)  exp  (by)  - i] 
- jp-  U - exp (-by) ] 

= -^  [i  - exp{-  (jp)y}]  . 


(3.31) 
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It  is  easy  to  observe  that  Eq.  (3.31)  turns  out  to  be 

identical  with  the  well-known  Janssen's  solution. 

(2)  Algebraic  solution.  From  Eq.  (3.28),  the 

lateral  pressure  pn  can  be  written 

n 


Substitute  the  constant  parameters  in  Eq.  (3.32)  to 
obtain  r 


The  series  included  in  the  square  brackets  is  a geometric 
series  with  a ratio  factor  smaller  than  1.  Summing 
the  geometric  series 


(3)  Comparison  with  Janssen's  solution.  Compare 
between  Eq.  (3.31)  and  Eq.  (3.35)  to  find  out  that  except 


i=l 


(3.32) 


pn  = Kpy [1  + (1-yK)  + ( 1-yK) 2+. . .+ ( 1-yK) n_1] . 


Pn  = Kpy  [ 


1- (1-yK) n 
1- (1-yK) 


= (1-yK)  n]  . 


(3.33) 


Using  Eq.  (3.27) 


(3.34) 


and  (3.33)  finally  becomes 


(3.35) 
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are  identical.  Now,  using  a limiting  process  where 
n-M»  (Ay+o)  , it  is  evident  that 

lim  r,  , 1.  yk.  -.n  r ,yk. 

n-H»  {1  + n{~—)Y}  = exp{-(H_)y}. 

(Ay->o) 

It  is  thus  shown  that  Eq.  (3.31)  and  Eq.  (3.35) 
become  practically  identical  for  sufficiently  small  Ay. 

iU 

Of  course,  it  was  shown  in  the  previous  section^Eq.  (3.31) 
is  identical  with  Janssen's  solution.  Thus  it  was  shown 
that  for  the  special  case  of  constant  parameters  both  the 
analytic  and  algebraic  models  become  identical  and  simply 
reduce  to  the  well-known  Janssen's  solution.  This  result 
is  very  important  as  it  defines  the  relationship  of  the 
new  theory  to  old  ones. 

3.  4.  2 Variable  Coefficients 

An  important  feature  of  the  new  theory  is  its 
great  flexibility  in  admitting  variable  parameters  of 
nearly  arbitrary  form.  A dual  test  of  the  models  for  a 
special  case  of  variable  coefficients  follows. 

(1)  Analytic  solution.  Explicit  integration  of 
the  analytic  solution  Eq.  (3.13)  is  very  sensitive  to  the 
nature  of  the  functions  y (y)  and  3 (y) , especially  the 
latter.  Unfortunately,  B(y)  is  bound  to  be  a very  com- 
plicated function  and  therefore  an  analytic  solution  is 
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seldom  practical.  However,  considerations  explained 
later  (see  section  4.3.1)  suggest  for  the  qualitative 
behavior  of  8(y),  under  certain  conditions,  the  nature 
of  a decreasing  hyperbolic  function.  The  simplest  form 
to  satisfy  this  requirement  is  a rational  function  of 
the  type 

; 6 (^  = IT^  (3‘36) 

where  n=— , A and  p are  constants.  The  ratio  -function 
P 

portion  of  3 (y)  is  then 

k(y>  = TlTcyT  (3'37) 

where  B and  C are  some  constants. 

The  density-function  is  an  increasing  function, 
and  for  simplicity  it  is  chosen  here  to  be  a linear 
function,  for  example 


y (y)  = Y0  + <$y  (3.38) 

where  yQ  and  6 are  constants.  Insert  Eqs.  (3 . 36) f (3. 38) 
into  Eq.  (3.13)  and  integrate  in  the  three  following  steps. 


(i)  Find 


(y)  dy  .J- 


nB 


(1+Cy) 


dy 


= Hi  Ind+Cy)  3/ 


3/  Constants  of  integration  may  be  deleted  throughout  this 
integration. 
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-j'y  (y ) exp  (Jt 


(ii)  Find  g(y)  =/y  (y)  exp  ( e (y)  dy ) dy=  / (y  +6y)  exp  [N  In  (1+Cy)  ] dy 


■h 


B 

where  N = v— ■ is  a constant.  Expanding  and  integrating  by 
parts,  one  finally  obtains 


N+l 


q(y)  = S-1+CY) [ (y  +<$y)  - M.^+Py.}.  ] 

gly;  C (N+l)  UYo  Y'  C (N+2 ) J* 


(iii)  Find  g (o)  = crN+'fy tYo  " CIMT1- 


Using  the  results  (i) , (ii)  and  (iii)  on  Eq.  (3.13), 
the  solution  is 


By 

p(y)=cWTJ 


o (1  1 , 6, N+l,  6 {l-(l+Cy)N+1}  , ,, 

U — + mjTt J • (3.39) 


(l+Cy)N+1  Yo  N+2  Yo  C(N+2) (l+Cy)N+1 


(2)  Algebraic  solution.  Use  the  same  3 and  y 
functions.  In  terms  of  k^,  Eq.  (3.28)  becomes 


n 


where 


p = k Ay  £ y ,,  .det  Sub.  [I  - rikAy] 
n 1 . , ' n+l-i  l 2 

i=l 


Yn+l-i  = Y0  + <5[n'  + 1 - i]Ay;  (n'=n  - i) 


X . . . 

n = - is  a constant 
P 


I is  an  n-by-n  identity  matrix 
B 


^n  1+Cn 1 Ay 


, and  k is  the  diagonal  matrix 


(3.40) 


36 


0 


0 


0 


• • 


0 


k = 


0 

0 


B 

1+C (n' -1) Ay 


B 

1+C (n ' -2) Ay 


0 


0 


0 


B 

1+C 1/2 Ay _ 


(3)  Comparison.  It  is  best  to  compare  between 
Eqs.  (3.39)  and  (3.40)  by  examining  their  graphs.  A 
numerical  example  was  worked  out,  and  the  results  are 
given  in  Fig. 

The  algebraic  solution  was  calculated  first  for 
large  increments  of  Ay  = 10*.  Twelve  discrete  points  of 
the  solution  were  obtained  and  joined  by  straight  lines. 

When  Ay  was  taken  smaller,  more  points  were  produced, 
and  the  algebraic  curve  tended  to  approach  the  analytic 
curve.  For  sufficiently  small  Ay  the  curves  practically 
coincided.  Thus,  it  is  shown  that  the  two  models  become 
identical,  in  the  limit,  for  variable  coefficients  as 
well. 

3.^.3  Relationship  of  the  Analytic  Model  to  the  Reimberts 1 

Semi-Empirical  Solution 

Further  investigation  in  this  connection  led  to  the 
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100  200 


Figure  3*1  Comparison  between  the  analytic  and  algebraic  solutions 
with  variable  coefficients. 
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following  important  result.  Starting  with  Eq.  (3.39), 
it  can  be  reduced  to  the  Reinberts ' semi-empirical 
solution  by  an  appropriate  choice  of  coefficients  and 
constants . 

First,  let  y (y)  = be  a constant  (i.e.  6=0), 
thereby  reducing  Eq.  (3.39)  to 

p(y),=  c(N+i)  [1  " (i+cy)- (N+1) ] . (3.41) 

Let  B = 2 kQ,  where  kQ  is  the  Rankine-Koenen-Caquot 
coefficient,  i.e. 


o l+sin<j> 


(3.42) 


where  <f>  is  the  angle  of  internal  friction  of  the  grain. 

Let  X = i y,  where  y is  the  coefficient  of  grain-wall 

yk 

friction,  and  let  C = — Then 


„ X B 2^  2ko 

N = p * c = — * Hk—  = 

(t£) 


Insert  the  values  of  N,  B,  and  C into  Eq.  (3.4  I)  which  thus 
reduces  to 


PYo  yko  -2 

P(y)  = ~ [l-{l+(_°)y}  2). 


(3.43) 
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It  is  now  easily  observed  that,  except  for  a 


difference  in  notation,  Eq.  (3.43)  is  precisely  identical 
to  the  Reimberts'  solution  reported  in  their  book  (43) 

(p.  35,  Eq.  (11))  in  the  form 

p = p [1-  (f  + l)-2] 

Mtiax  A 


where 

z 

depth  below  the  edge  of  the  vertical  wall 

Pz  - 

lateral  pressure  exerted  on  the  wall  by 

p 

rmax 

the  grain  at  depth  z 

6D  . 

4tg<j>'  ' where 

6 

density  of  grain 

D 

interior  diameter  of  cylindrical  bin 

r 

angle  of  friction  of  grain  on  wall 

A 

D h . . . 

— T is  a constant 

4tg+-tg2(J-  - f)  3 

h 

height  of  a conical  surcharge  (zero  in  this 
comparison) 

4> 

angle  of  internal  friction  of  grain 

+ 2,TT  = 1"sin(l)  = 

tg  4 2 l+sin(j)  o 
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Discussion.  The  Reimberts ' solution  was  derived 


on  the  basis  of  both  theoretical  considerations  and 
empirical  data  gathered  over  years  of  experimental 
investigations  of  dynamic  effects  in  grain  bins.  The 
derivation  of  their  method  is  different  from  the  derivation 
of  the  theories  of  this  study,  which  is  strictly  theoretical. 
However,  it  is  most  interesting  to  note  that  our  analysis 
supplies  a fundamental  answer  to  the  understanding  of  the 
Reimberts*  method  and  results,  not  known  before. 

Assuming  that  the  Reimberts*  solution  is  in  good 
agreement  with  their  experimental  findings,  the  result  of 
Eq.  (3.43)  determines  the  exact  behavior  of  the  ratio- 
function  and  wall-friction  parameter  during  their  experiments 
The  essence  of  the  Reimberts*  method  in  relation  to  ours 
can  be  summarized  as  follows: 

(i)  the  density-function  yQ  is  a constant 

(ii)  the  wall-friction  parameter  ^ y is  a constant 

(iii)  the  ratio-function  is  a decreasing  hyperbolic 

rational  function  of  the  form 

2k 


k(y) 


1+  (■ 


yk 


R 


)y 


(iv)  the  lateral  pressure  is  calculated  according  to 
our  general  analytic  solution  in  one-dimension, 

Eq.  (3.13),  using  the  above  factors  as  coefficients 
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These  conclusions  are  very  significant.  First, 
they  establish  explicit  mathematical  expressions  for  a 
dynamic  ratio- function  and  a dynamic  wall  parameter  X. 

It  is  confirmed  mathematically  for  the  first  time  that  the 
dynamic  ratio-function  is  a hyperbolically  decreasing 
function.  Also,  that  under  dynamic  conditions  the 
friction-parameter  is  one-half  the  static  friction  co- 
efficient. It  should  be  noted  that  without  the  generalized 
solution  it  would  have  been  next  to  impossible  to  gain 
insight  into  these  complex  interrelationships  under 
dynamic  conditions.  Furthermore,  at  this  point,  it  was 
established  that  both  the  Janssen's  and  the  Reimberts' 
solutions  --  thought  until  now  to  be  independent  --  are 
not  independent,  but  are  simply  special  cases  of  the 
same  prototypical  model,  Eq.  (3.13).  Many  more  types  of 
solutions,  static  and  dynamic,  can  be  derived  from  our 
model  as  shown  in  Table  3.1.  Through  the  use  of  the  digital 
computer,  these  models  can  be  realized  and  then  exploited 
for  both  research  and  design  purposes. 
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TABLE  3.1.  Analytic  solution  as  a prototype  of  other  solutions 
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CHAPTER  4 


THE  CHARACTERISTIC  FUNCTIONS 
4.1  INTRODUCTION 


The 

characteristic  functions 

are : 

(1) 

The 

density- function 

Y(y) 

(2) 

The 

ratio- function 

My) 

(3) 

The 

.friction- function 

X (y)  • 

They  determine  the  3 and  y functions  and  as  such 
are  the  building  blocks  of  the  basic  model.  They  have  to  be 
explicitly  specified  before  a model  is  fully  characterized 
and  implemented.  As  mentioned  before,  our  models  are  very 
general  and  will  admit  arbitrary  characteristic  functions. 

In  other  words,  there  is  absolutely  no  theoretical  limitations 
on  these  functions  and  various  models  can  be  produced  by 
arbitrary  choices.  Of  course,  some  choices  may  lead  to  mathe- 
matical constructs  that  have  no  physical  counterparts.  These 
are  usually  discarded.  However,  many  do  lead  to  interesting 
constructs  that  may  be  considered  true  models.  The  ultimate 
ease  of  realizing  the  models  that  develops  from  their  imple- 
mentation on  the  digital  computer,  makes  it  reasonable  to 
experiment  with  a variety  of  models  based  on  theoretically 
prescribed  characteristic  functions.  This  process  of  experi- 
mentation is  called  digital  simulation  and  will  be  discussed 
later.  Ultimately,  empirical  characteristic  functions  have  to 
be  fed  into  the  basic  model  to  establish  practical  solutions. 
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A discussion  of  these  functions  under  static  and 
dynamic  conditions  follows. 

4 . 2 THE  DENSITY-FUNCTION 


The  apparent  density  of  grain  stored  in  bins  is  defined 

* as  the  weight  of  grain  bounded  by  a unit  volume.  Most  types 

of  grain  have  rather  definite  density  values,  although  the 
r 

density  value  of  a particular  type  of  grain  may  vary  within 
a range  of  some  ±10  percent. 

Variations  in  grain  density  may  be  due  to  factors  such 
as:  moisture  content,  compactness  (looseness),  pre-sure, 
vibrations,  settlements  and  so  on.  Under  normal  storage 
conditions  the  density  value  is  bounded  between  some  minimum 
and  maximum  values  which  may  differ  by  10  to  20  percent. 

If  not  considered  merely  as  a constant,  the  one- 
dimensional density-function  is  expected  to  be: 

(1)  continuous  and  smooth 

(2)  monotonically  increasing  with  depth 

(3)  strictly  positive 

(4)  bounded  between  y . = y and  y , empirical  values. 

min  o max 

The  simplest  function  to  satisfy  these  requirements  is 
a linear  function 


y0  + 


y(y)  = 


.Y 


max 


for  0 < y < 


for  Y > Hd 


(4.1) 
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where  6 is  a constant  defined  by 


Y “Y 
* 'max  'o 

0 = 75 . 

Hd 

If  H is  considerably  larger  than  H^,  it  is  preferable 
to  choose  a nonlinear  function.  It  seems  reasonable  to  expect 
that  the  increase  in  density  versus  depth  is  mainly  proportional 
to  the  vertical  pressure.  Therefore  Y (y)  is  assumed  to 
increase  exponentially  according  to  Janssen-type  (static) 
pressure,  i.e. 

yk 

Y(y)  = Yq  + (Ymax  - Yq)  [ 1-exp {-(-jp)y}] . (4.2) 

For  most  practical  cases  in  grain  bins  it  is  sufficient 

to  assume  that  the  density-function  is  a constant,  y , 

ave 

or  a linear  function  such  as  Eq. (4.1). 

Eq.  (4.2)  is  particularly  well  fitted  to  the  density- 
function  in  silage  silos. 

The  density-function  is  independent  of  dynamic  conditions. 


4 . 3 THE  RATIO-FUNCTION 


The  one-dimensional  ratio-function  k(y)  is  defined  as 
the  ratio  of  the  lateral  pressure  to  the  vertical  pressure 
at  any  depth  y,  i.e. 


(4.3) 


46 


k(y)  is  essentially  an  empirical  factor.  In  order 
to  investigate  the  nature  of  k(y)  it  is  necessary  to  rely  on 
experimental  findings  reported  from  various  sources.  It 
seems  fair  to  summarize  these  findings  as  follows: 

(1)  Under  static  conditions  the  actual  lateral  pressure 
L closely  agree  with  Janssen's  curves. 

(2)  Under  dynamic  conditions  the  actual  lateral  pressure 
may  exceed  Janssen's  curves  (by  tens  to  hundreds  percents) 

until  they  achieve  a maximum  at  some  depth  above  the  bottom, 
and  then  decrease  rapidly.  The  behavior  is  experimentally 
not  clear  in  this  lower  region. 


4.3.1  Static  Conditions 


Under  static  conditions  no  grain  is  added  or  removed 
from  the  system,  and  it  is  assumed  that  all  parameters  are 
time- independent  and  that  the  walls  do  not  yield  appreciably 
under  the  action  of  the  lateral  pressure. 

Janssen  (1895)  who  first  introduced  the  ratio-f actor , 
assumed  that  it  is  an  empirical  constant.  Koenen  (1896) 
suggested  to  use  the  theoretical  "coefficient  of  earth 
pressure  at  rest"  to  represent  this  constant.  Thus  Rankine's 
coefficient  was  borrowed  for  use,  i.e. 


k 


o 


l-sin(f) 

1+sincf) 


(4.4) 
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where  <f>  is  the  angle  of  internal  friction  of  the  grain, 
treated  as  an  empirical  constant.  The  above  form  of  kQ 
was  adopted  by  most  methods  and  is  found  to  be  sufficient 
for  the  case  of  static  pressure.  It  is  therefore  proposed 
to  adopt  it  in  the  present  analysis  also.  However,  there 
arise  some  important  notes. 

By  Eq.  (4.4),  kQ  is  a function  of  <f> , the  angle  of  internal 
friction  of  grain  which  is  defined  as  "the  angle  whose  tangent 
equals  the  ratio  between  the  shearing  resistance  per  unit  area 
to  the  corresponding  normal  stress  in  non-cohesive  grain. " 

4)  is  actually  a measure  of  the  stability  of  the  grain.  It  is 
frequently  confused  with  the  angle  of  repose  6 which  indeed 
often  has  the  same  value  and  also  serves  as  a measure  of 
stability . 

4>  is  only  approximately  equal  to  0.  For  most  materials, 

although  not  all,  it  is  slightly  larger  than  0.  Now,  since 

4>  depends  on  the  local  normal  stress  (pressure)  , it  is  not 

constant  and  therefore  kQ  is  not  constant.  Whenever  4)  is 

considered  constant  the  effect  on  k is  similar  to  using  0. 

o 

This  practice  is  however  always  on  the  safe  side,  since  the 
resulting  kQ  is  larger  than  it  would  have  come  out  otherwise, 
namely 


k 


o max 


l-sin0 

l+sin0 


(4.5) 
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The  influence  of  the  variable  <J>  (y)  on  kQ  will  now 
be  investigated.  First,  establish  an  approximation  to  the 
function  (y)  . <f>  is  primarily  dependent  on  the  density 

(and  the  vertical  pressure) . If  the  relationship  is  assumed 
to  be  proportional  then,  qualitatively,  <J)  (y)  behaves  like 
Y (y) r as  was  established  experimentally  by  Platonow  and 
Kowtun  (1959).  Using  Eq.(4.2),  <f>  (y)  can  be  written 

<f>(y)  = 4>  + (4>  ) [l-exp{-( — °.Tna-)y}]  (4.6) 

o max  o p 

where 

<f>Q  - angle  of  repose,  0 

<Pjnax=  a(PQ'  where  a is  an  empirical  coefficient  depending, 
among  other  factors,  on  H;  i.e.,  a=2  approximately 
in  deep  grain  bins,  according  to  the  findings  of 
Platonow  and  Kowtun  (1959). 

k is  defined  by  Eq.(4.5). 

o max  -1 

Then,  substituting  Eq.  (4 . 6)  in  Eq.(4.4),  ko(y)  is 
finally  obtained  in  the  form 


ko(y) 


1-sin  [0 (2-exp{- 


1+sin  [0  (2-exp{- 


y (l-sin0) 
p (l+sin0 ) 
y (l-sin0 ) 
P (1+sin® 


y})] 

y})] 


(4.7) 


and  kQ(y)  is  thereby  expressible  as  a function  of  y, 
involving  0,  y and  p as  parameters.  The  behavior  of  the 
function  ^Q(y)  may  be  seen  by  using  the  trigonometric 
identity 
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(4.8) 


l-sin<f) 

1+sincj) 


. 2 ,TT 

= tan  (j  - 


Let  (J  ”~)  = ip  and  examine  \p . By  Eq . (4.6),  ^ is  increasing 
versus  y,  and  therefore  is  decreasing.  The  tangent  of  the 
decreasing  argument  \p , where  0 c ip  < is  decreasing  versus 
y and  its  square  is  even  more  so.  It  follows  that  kc (y)  is 
a decreasing  function  of  y,  convex  toward  the  y-axis.  Such 
behavior  can  well  be  represented  by  the  appropriate  choice 
of  a hyperbola  such  as  the  type  used  in  section  3.3.2. 

4.3.2  Dynamic  Conditions 

It  is  generally  accepted  that  dynamic  effects  occur 
during  charging  and  especially  discharging.  The  discussion 
below  is  confined  to  discharging  alone. 

As  mentioned  before,  the  determination  of  k(y) 
should  be  done  experimentally.  In  fact,  under  dynamic  con- 
ditions the  ratio-function  must  be  time-dependent,  k(y,t), 
and  this  further  complicates  the  picture.  In  the  absence 
of  real  experimental  dynamic  data  the  discussion  below  is 
somewhat  speculative. 

Suppose  that  as  a result  of  a discharging  operation, 
the  walls  deflect  to  the  extent  that  the  contact  between 
the  grain  and  the  walls  is  reduced  to  a minimum.  At  this 
particular  instant  two  main  phenomena  may  occur: 
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(1)  Since  the  friction  effect  approaches  zero,  the 
vertical  pressure,  and  therefore  the  lateral  pressure  as  well, 
rapidly  increase.  If  the  deflection  of  the  walls  is  sufficiently 
large,  so  that  the  system  can  be  looked  upon  as  if  the  walls 
were  completely  removed  (for  a moment) , the  column  of  grain 
tends  to  collapse  and  transform  into  the  characteristic 

pile  (see  the  discussion  in  the  section  on  topological  models) . 

(2)  At  the  same  instant  the  walls,  being  released 
from  pressure,  tend  to  return.  The  assumption  is  that  the 
maximum  lateral  pressure  develops  when  the  returning  walls 
and  the  nearly  collapsing  grain  interact. 

"Dynamic  conditions"  will  be  regarded  as  conditions 
under  which  the  two  following  requirements  are  satisfied: 

(1)  The  lateral  pressure  is,  qualitatively  at  least, 
in  agreement  with  the  experimental  findings  described  in 

the  previous  section. 

(2)  The  walls  deflect  appreciably. 

Further  studies,  both  digital  simulations  and  experi- 
mentation, should  be  made  to  establish  reliable  k(y,t) 

Is  functions  under  various  dynamic  conditions.  In  Chapter  7 
a theoretical  attempt  is  made  to  characterize  a dynamic 

I 4 

ratio- function. 

' 
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4.4 


THE  FRICTION-FUNCTION 


The  friction-function  X (y)  is  defined  as  the  ratio 
of  the  frictional  stress  at  the  wall  to  the  lateral  pressure, 
at  any  depth,  namely 

X(y)=|gf.  (4.9) 

The  coefficient  of  friction  of  the  grain  on  the  wall 
y is  basically  dependent  on  the  type  of  the  grain,  the 
material  of  the  wall,  and  the  degree  of  smoothness  (roughness) 
of  the  wall  surface.  Other  factors  such  as:  variations  of 
the  moisture  content  of  the  grain,  and  the  motion  of  the 
grain,  may  affect  the  coefficient  of  friction.  Nevertheless, 
this  factor  is  normally  considered  as  a constant. 

The  wall-friction  parameter  X was  introduced  in  order 
to  accommodate  variations  in  the  friction  effect,  either  as 
a function  of  the  depth  coordinate  y,  or  independent  of  y, 
but  as  a variable  parameter  depending  on  other  factors. 
Therefore , 

X=  vy  (4.10) 

where  v is  either  (i)  a constant,  (ii)  a function  of  y,  or 
(iii)  a variable,  but  independent  of  y.  The  flexibility  in 
the  definition  of  v allows  to  take  into  account  many  possible 
variations  in  the  friction  effect.  v is  made  to  be  bounded 
between  0 and  1,  such  that  X is  bounded  between  0 and  y. 
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4.4.1  Static  Conditions 


The  friction  effect  attains  its  maximum  influence 
under  static  conditions.  Therefore  let  v=l  such  that 


A 


A 


max 


V- 


(4.11) 


4.4.2  Dynamic  Conditions 


Since  v is  bounded  between  0 and  1,  it  may  prove 
useful  to  express  it  in  the  form  of  a trigo-ometric  function, 
i.e.  the  sine  function.  Let 

v(y)  = sin  cy  (4.12) 


where  c is  a parameter. 

This  representation  may  be  appropriate  especially 
for  a case  where  the  wall  vibrates  under  dynamic  conditions. 
When  the  wall  vibrates,  it  is  likely  to  deflect  in  the  form 
of  a sine  curve.  The  mode  of  the  vibrations  is  not  known 
(involved  in  c) , but  can  be  estimated  for  simple  cases.  It 
is  assumed  that  the  wall  is  much  more  flexible  than  the  mass 
s of  the  grain  is,  and  that  the  grain  mass  cannot  continuously 
take  on  the  form  of  the  deflecting  wall.  Therefore  the  grain 

* 

is  in  very  loose  contact  with  the  wall,  or  in  no  contact 
at  all,  along  the  intervals  where  the  wall  is  convex  outwards. 
Therefore  define  v(y)  as  follows: 
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sin  cy 


for  sin  cy  > 0 


v (y) 


(4.13) 


0 


for  sin  cy  < 0 


From  the  nature  of  the  sine  curve  and  definition  (4.13), 
it  is  evident  that  regardless  of  the  mode  of  the  vibrations 
v(y)  vanishes  intermittently  over  a cumulative  length  of 
one-half  of  the  total  depth  of  the  bin. 


factor,  it  is  not  necessary  to  find  the  exact  mode  of  the 
deflections.  The  mode  can  be  taken  for  convenience  as  the 
arbitrary  interval  Ay,  according  to  which  the  numerical 
calculations  are  carried  out.  This  choice  produces  a v(y) 
that  vanishes  at  every  other  interval  and  varies  between  0 
and  1 according  to  a sine  curve  at  every  other  interval.  For 
numerical  calculations  of  sufficiently  small  Ay,  v (y)  ap- 
proaches the  form  of  a step-function.  Therefore  it  can  be 
written  as  follows. 


Since  the  friction  effect  is  important  as  a cumulative 


sin  cy  for  sin  cy  > 0 


v (y) 


(4.14) 


for  sin  cy  < 0 


where 


ave.  sin  cy  = — 

2 77 


2 
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Therefore,  under  "simple"  dynamic  conditions  and  a choice 
of  sufficiently  small  increment,  v (y)  may  be  chosen  to  a 
step-function  of  the  form 


v (y)=< 


for  i odd  integer 


for  i even  integer 


(4.15) 


where  i is  related  to  Eq.  (3.28),  i.e.  i=y/Ay.  Accordingly 
the  friction-function  can  be  expressed  as  a step-function 
of  the  form 


for  i odd  integer 


for  i even  integer  . 


(4.16) 


4.4.3  Practical  Recommendations 


Eq.  (4.16)  may  suggest  to  choose  a constant  friction- 
function  which  is  the  average  of  the  intermittent  values  0 
and  , namely 


> 


The  choice  of  constants  smaller  than  this  is  justified 
for  cases  where  the  assumption  of  simple  sinusidal  deflec- 
tion is  not  valid.  In  any  case  of  uncertainty  it  is  ad- 
visable to  choose  smaller  values  for  the  parameter  v,  in- 
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eluding  the  possibility  of  a zero  value.  When  v=0,  the 
friction  effect  is  completely  neglected  and  thus  the  lateral 
and  bottom  pressures  attain  maximum  values.  However,  to 
calculate  the  vertical  load  on  the  walls,  X =y  should 
be  used. 

In  summary,  there  may  exist  three  cases  of  significant 
difference  (among  many  possible  intermediate  cases) : 

(1)  Static  conditions.  v is  constant  = 1.  The 
friction  parameter  attains  its  upper  bound,  ^max=y« 

(2)  "Simple"  dynamic  conditions.  v(y)  is  a step- 
function  of  the  form  given  by  Eq.  (4.15);  or,  more  prac- 
tically, v can  be  taken  as  an  averaged  constant,  namely, 
v = | and  therefore  X = -jy . 


(3) 


Complex  dynamic  conditions.  y = 0 and  X = 0. 


For  design  purposes  the  following  is  suggested: 

(1)  For  rigid  structures  with  special  discharging 
devices  (such  as  a central  perforated  pipe)  use 

X = X = y . 
max 

(2)  For  rigid  structures  without  special  discharging 
devices,  or  for  semi-rigid  structures  with  such  devices,  use 


X = 3 v' 

(3)  For  non-rigid  structures  with  hazardous  possi- 
bilities for  dynamic  effects,  use 
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A-*0,  or  A=0. 


Notes . The  following  are  general  examples  of 
different  types  of  grain  bins  that  may  occur  in  practice, 
with  respect  to  the  degree  of  "rigidity"  as  used  in  the 
above  recommendations.  For  more  accurate  determination  of 
any  specific  design  case,  a complete  structural  analysis 
should  be  made. 

Rigid  structures  are  structures  such  as  clusters 
of  monolithic  reinforced  concrete  bins  with  rigid  connections 
to  a rigid  slab  foundation,  based  on  a high-quality  soil 
and  with  rigid  connections  at  the  tops. 

Semi-rigid  structures  are  single,  exceptionally 
tall  reinforced  concrete  bins;  clusters  of  steel  bins. 

Non-rigid  structures  are  single,  tall,  steel  bins 
with  a non-rigid  connection  to  the  foundation,  based  on 
an  inferior  soil  and  subjected  to  frequent  gusts. 

Proper  discharging  devices  are  such  as  those 
reported  in  the  literature  and  tested  in  operation  in  many 
countries.  However,  the  design  and  installation  of  such 
devices,  do  not  necessarily  guarantee  that  the  dynamic 
effects  will  be  eliminated.  Therefore,  a conservative 
choice  of  the  values  of  the  parameter  v is  encouraged. 


CHAPTER  5 


NUMERICAL  AND  COMPUTER  MODELS 

5.1  NUMERICAL  MODELS  BASED  ON  THE  ANALYTIC  MODEL 

The  numerical  models  of  this  category  are  based  on 
the  procedure  described  in  section  3.2.  The  model  is 
written  first  in  the  form  of  Eqs.  (3.18)  and  (3.19)f  B (y) 
being  specified  by  the  particular  k(y)  and  X (y)  chosen 
for  the  model.  Then  Eq.  (3.18)  is  programmed  into  the 
FUNCTION  subprogram  named  FUN.  The  control  program  calls 
the  fourth-order  Runge-Kutta  integration  subroutine  RK2 
(identical  with  subroutine  RUNKUT  given  in  section  3.2). 
Subroutine  RK2  returns  the  discrete  solution  for  the  vertical 
pressure  q(y).  The  control  program  calculates  the  lateral 
pressure  p(y)  by  multiplying  q(y)  by  k(y). 

5.1.1  Examples  and  Results 

(1)  Example  No.  1.  The  characteristic  functions 

are: 

Density- function: 

y = 45.0  lb/cu. ft.  (constant). 

Ratio- function: 

2k 

k(y)  = — 

1+0 . 06y 
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where , 

k = 0.333. 
o 

Friction- function: 

X = y = 0.450  (constant). 
The  system  solved  is: 


, 2k 

dq  _j_  y o 

dy  Y “ R • 1+0. 06y  - q 

q (0 ) = 0 

Geometry: 

H = 100.00  ft. 

R = 4.99  ft. 

Ay  = 1.00  ft. 

The  final  results  give  the  lateral  pressure,  p(y), 

2 

in  tabular  form  in  both  kg/m  and  Ib/sq.ft.  The  same  results 
are  also  plotted  at  2.0  ft.  and  1.0  ft.  intervals. 
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NUMERICAL  SOLUTION  OF  THE 
ANALYTIC  MODEL 

Example  No.  1 


A] 


1) 


DIMENSIO 
EXTERNAL 

1 FORMAT!3 

2 FORMAT (1 
11H  tlOX  ,2 
21H, 12X,  1 

3 FORMATd 
10  RE  AO  (5, 

2D  WR  IT  E { 6, 

CALL  RK2 
STEP  =F  LO 
X=XO 

DO  30  1 = 
XX  ( I ) = X 
X=  X+  STEP 
A!  I ) = A(I  ) 
B=A ( I)*4 
NR  ITE  ( 6 , 
751  F3RMAT(  1 
WR I TE ( 6* 
CALL  AMA 
CALL  AMI 
CALL  PL 
CALL  PLO 
VA  L M A X = 2 
CALL  PLO 
CALL  PLO 
WRITE! 6, 
75b  FORMAT! • 
4P  CONTINUE 
RETURN 
END 


A ( 500  ) ,X|X(  500) 

FUN 

F10. 0,215) 

HI, 7X,44HS 
HH= , F7  .3, 2 
HX  , 18  X , 4H  Y 
H,  1 OX  ,F5 . 2 
. XO,YO,H, 

2)  H, XC, YO 

!FUN,H,  X0,|Y0,  JNT,  I EN|T  , A ) 
AT! JNT) *H 


OLUTION  OF 
X,  3H X 0 = , F 7 
IX)//) 

, LOX  ,E  1 5 
JNT,  I ENT 


1,1  ENT 


! !0.666) 

.88  2 

751)  X ,B 
0X,F10. 5,1 
3)  X , A ! I ) 
X!  A, IENT  ,V 
N!  A, IENT, V 
{ 1 , 1 ENT, I 
IT!0,IENT,1 
.0*VALMAX 
T ! 1 , I ENT  , 1 
T ( 0 , IENT  , 1 
750)  IENT 
TOTAL  NO. 


OT 


DY/DX=FUN( 
L 3 ,2  X , 3 H YO 


8) 


/(I .0+0.06 


*XX! I) ) ) 


0X,F15.6) 

ALMAX , I SUB ) 
ALM IN, I SUB ) 
, VALMAX , VA 
, VALMAX, VA 


VAL 


, VA  LMA  X , 

, VALMAX, VA 

OF  INCREM 


LM  IN, XX, A ) 
LM  IN  ,XX  , A ) 


MIN, XX, A) 
LMIN, XX, A) 

ENT S = ' , I 5 ) 


X ,Y)  BY  R 
, F7  .3// 


<2  SUBROUTINE// 


FUNCTION 

F UN ( X , Y ) 

FUN=45 • 0 
RE  TURN 
END 

-(0.450/4. 

90)*( (0.66 

6 ) / ( 1 . + ( 0 . 

06* X ) ) ) * Y 
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SOL 


UTION  OF 
H=  1.0C0 


DY/0X=FUN( 

X,Y)  BY 

RK 

2 SUBR0UTI 

X 

o 

II 

o 

a 

o 

Y0  = 

0 

.0 

Lateral  pressure  p(y) 

Depth  (ft.)  (kg./m.2  and  lb./sq.  ft.)  Depth  (ft.) 


i.COOOO 

142.  1 

68350 

26.00000 

1.0C 

0.2 

9120926E  0 

2 

26.00 

2. CC0C0 

261.2 

59766 

27.00000 

2.00 

0.5 

3514923E  0 

2 

0?  7?  99  6 7 8 9 

3.00000 

361.9 

90967 

28. COOOO 

3.00 

0.7 

4148117E  0 

2 

28.  CC 

4. C0000 

447.  9 

3 53  03 

29.00000 

4.CC 

0.9 

1752426E  0 

2 

o29  jOO  6 7 8 9 

5.CCCC0 

521.8 

37158 

30.00000 

5.0C 

0.1 

0689005E  0 

3 

30.00 

6.00000 

585.8 

34961 

31  ..COOOO 

6.00 

0.1 

19999C1E  0 

3 

31. OC 

7.00000 

641.6 

123  05 

32.C00C0 

7.00 

0.1 

314241 CE  C 

3 

32.00 

8.C000O 

690.5 

11963 

33.00000 

8.CC 

0.1 

4 144043E  0 

3 

33.00 

9.00000 

733.6 

12793 

34. COOOO 

9.00 

0.1 

5026892E  0 

3 

34.00 

10.00000 

771.7 

91016 

35. CC000 

10.00 

0.1 

5 80891 3E  0 

3 

35. OC 

11.C0000 

805.  7 

64404 

36.00000 

11. oc 

0.1 

650480 3E  0 

3 

36  .00 

12.CC0C0 

836.1 

23291 

37. COOOO 

12.00 

0.1 

7126656E  C 

3 

37.00 

13.00000 

863.3 

59619 

38. COOOO 

13.00 

0.  1 

76  8455  CE  0 

3 

38. OC 

14.C0000 

887.  8 

83  545 

39.00000 

14. OC 

0.1 

8 186885E  0 

3 

39.00 

15.CC0C0 

9 10.0 

41504 

40.00000 

15. OC 

0.  1 

8640755E  0 

3 

40.00 

16.00000 

930.1 

25000 

41.00000 

16.00 

C.l 

9C52132E  0 

3 

41. OC 

17.00000 

948.3 

83789 

42.00000 

17.00 

0.1 

9426135E  C 

3 

42.00 

18.C00C0 

965.0 

29297 

43.00000 

18.00 

0.1 

9767093 E 0 

3 

43.00 

19.00000 

980.2 

^536  1 

44. COOCO 

19.00 

C.  2 

0078767E  0 

3 

44.00 

20.00000 

994.1 

89209 

45. COOCO 

20.00 

0.2 

0364384E  0 

3 

45.00 

21.C0000 

10C6. 9 

9707C 

46.00000 

21. OC 

0.2 

0626735E  0 

3 

46.00 

22.C00C0 

1018.7 

88818 

47.00000 

22.00 

0.2 

0868268E  0 

3 

47.00 

23.00000 

1029.6 

66504 

48. COOOO 

23.00 

0.2 

1C91C81E  0 

3 

48. OC 

24.00000 

1039. 7 

22412 

49.00000 

24.00 

0.2 

1297060E  0 

3 

49.00 

25. C0C00 

1049.0 

35645 

50.00000 

25. OC 

0.2 

1487827E  0|3 

50.00 

Lateral  pressure  pCy) 
(kg./m.2  and  lb./sq.  ft.) 

1057.6177002 
0.21 664 833 E 03 

1065.708984 

34  2 34  5 - 

1073. 1 86279 
0.21982516E  03 

1080.159180 

345  60^2212534  51  8 00  1 214S: 

1086. 670166 
0.2225871 3E  03 

1092. 760010 
0. 2238  3450E  03 

1098.462402 
0.22500256E  03 

1103.810303 
0. 2 2609798E  03 

1108. 830566 
0.22712634E  03 

1113.550537 
0. 2280931 2E  03 

1117.992187 
0.229CC293E  03 

1122.177246 
0 .2298602  CE  03 

1126.125244 
0.23066884E  03 

1129.852539 
0.23143231 E 03 

1133.375244 
0.2321 5392E  03 

1136. 7C8008 
0.23283655E  03 

1139.864258 
0.23343306E  03 

1142.854980 
0 .2 34 05567E  03 

1145. 692383 
0.23467688E  03 

1148.386475 
0. 23522  87CE  03 

1 1 50 • 9 V6045 
0.23575299E  03 

11 53.380371 
0.23625166E  03 

1155. 696533 
0. 23672606E  03 

1157.902588 
0.23717796E  03 

1160.005127 
0.2376C86CE  0|3 
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Lateral  pressure  p(y) 

Depth  (ft.)  (kg./m.?  and  lb./sq.  ft,)  Depth  (ft.) 


51. CC000 
51.  OC 

52.  COCCO 

52.00 

53.00000 
53  .00 

54.00000 

54.00 

55. COOCO 

55.00 

56.00000 

56.00 

57.00000 

57.00 

58. 00000 
58. OC 

59.00000 

59.00 

60.00000 

0^?  Q ^ Q 9 6 7 8 9 01 

61.00000 

61.00 
62.C0000 


1162. 

0 • 2 380 19  30F 
1163. 

0.23841 124E 
1165. 

0. 2387854 5E 

1167. 

0 .23914323E 

1169. 

0. 23948524E 

1170. 

0.2398123 8E 

1172. 

0 .240 1 256  IE 

1173. 

0. 24042554E 

1175. 

0.24C71300E 

1176. 

2 3 4 5 8 

1177. 

0.24 12531 6E 
1179. 


>10010 

i3 

>23584 

>3 

’50732 

>3 

> 97070 
) 3 

66748 

03 

,763916 

03 

2 93213 

03 

757324 

03 

150645 

03 

5 06836 

Qp. 
797852 
0 3 

036621 


12  3 4 5 6 


O62l00  4 7 8 9 

01234560^2 

41 5 0691 E sO 

3 1 2 3 4 

63.00000 

1180.2 

26563 

63.00 

0.2 

4175C67E  0 

3 

64. COOOO 

11 81.3 

69873 

64.00 

0.2 

4198486E  0 

3 

65. COOOO 

1182.4 

68994 

65.00 

0.2 

4220999E  0 

3 

66.00000 

1183  .5 

25635 

66.00 

0.2 

4242641 E 0 

3 

67.00000 

1184. 5 

42725 

67.00 

0.2 

4263475E  0 

3 

68. COOOO 

1185. 5 

21240 

68.  OC 

0.2 

4283516E  0 

3 

69.00000 

1186.4 

63623 

69.00 

0.2 

43  02  82  OE  0 

3 

70.00000 

1187.3 

71582 

70.00 

0.2 

4 2 2 1 42  2E  0 

3 

71.00000 

11 88.2 

46338 

71.00 

0.2 

43  39340E  0 

3 

72. COOOO 

1189.0 

89844 

72.00 

0.2 

435661 3E  0 

3 

73.00000 

1189.9 

03076 

73  .00 

0.2 

43  73274E  0 

3 

74. COOOO 

11 90. 6 

88232 

74.00 

0.2 

43  8935  5E  0 

3 

75. COOCO 

1191.4 

45557 

75.00 

0.2 

4404869E  C 

3 

76.00000 

76.00 

77.00000 

77.00 

78. COOOO 

78.00 

79.00000 

7 9.00 

80.00000 

80.00 
81.00000 

81.00 

82.00000 

8 2.00 

83.00000 
83  .00 

84.00000 

84.00 

85. COOCO 
8 5.00 

86.00000 

86.00 

87.00000 

87.00 

88. COOOO 

88.00 

89.00000 

89.00 

90.00000 

90.00 

91.00000 

91.00 

92. COOCO 

92.00 

93.00000 

9 0??  ?99  6 7 8 9 I 0 I 2 

94.  COOOO 

94.00 

95.  COOCO 

9 o95i0G  6 7 8 9 0 12 

96.00000 

96.00 

97. 00000 

97.00 

98. COOOO 
98. OC 

99.00000 

99.00 

l 00. COOOO 
$ $ $ $ $ 


Lateral  pressure  p(y) 
(kg./m.2  and  lb./sq.  ft.) 

1192.1(77002 
0.2441 9852E  03 

1192. 883057 
0.2 443431 2E  03 

1193.565674 
0.24448296E  03 

1194.225342 
0. 24461809E  03 

1194. 863037 
0 • 244  74  8 72E  03 

1 195. 4|8C469 
0 • 2448751 5E  03 

1196.077148 
0.2  449974 1 E 03 

1196.654297 
0.24511563E  03 

1197.213867 
0.2452  30  2 2E  03 

1197. 755371 
0.24534117E  03 

1198.279297 
0. 24544  847E  03 

1198. 787842 
0.245552  64E  03 

1 199.279785 
0.24565 34 3E  03 

1199.757324 
0.245 75 122E  03 

1200.220459 
0.245  846 C8E  03 

1200.669922 
0.24593814E  03 

1201. 105225 
0.24602733E  03 

1201 .528076 

34  2 34  56 

1201. 938477 
0 .246  1 9798E  03 

1202. 3371  58 

3 4 5 6©  ^ 2|4  62  796  6{f  8 0|3  1 2 3 4 5 6 

1202.723877 
0. 24635889E  03 

1203.099854 
0. 2464359CE  03 

1203.  464  844 
0.24651067E  03 

1203 .8  19824 
0.246583 3 7E  03 

1204. 1 64795 
0.2466  5404E  03 
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PLOTTING 
DEPENDENT 
Depth  ( 

SY  MB  CL 
VARIABL 
ft.) 

I NDE  PENT 

ENT 

VARIABL 

E 

0 .0 

2.000000 

OE 

00 

4. CCOOOO 

OE 

00 

6 • OC  0000 

CE 

00 

8.000000 

OE 

OC 

1.000000 

OE 

01 

i . 200000 

OE 

01 

1.4CC00C 

OE 

01 

1.600000 

OE 

01 

1.80000C 

OE 

Cl 

2.000000 

OE 

01 

2.200000 

OE 

01 

2.4C0C0C 

CE 

01 

2.60C00C 

CE 

01 

2.80000C 

OE 

Cl 

3.000000 

OE 

01 

3.200000 

OE 

01 

3.4C0000 

OE 

01 

3.60000C 

CE 

01 

3.800000 

OE 

01 

4.000000 

OE 

01 

4. 200000 

OE 

01 

4.4CC000 

OE 

01 

4.60000C 

OE 

01 

4.800000 

OE 

01 

5.000000 

OE 

01 

5.2CC000 

OE 

01 

5.4CC00C 

OE 

01 

5.600000 

OE 

Cl 

5.800000 

OE 

Cl 

6.000000 

OE 

01 

6. 2 COO 00 

OE 

01 

6.4C0000 

OE 

01 

6.60000C 

OE 

01 

6.800000 

OE 

01 

7.000000 

OE 

01 

7.2CC000 

OE 

0 1 

7.4CC000 

OE 

01 

7.6CC000 

OE 

01 

7.800000 

OE 

01 

8.000000 

OE 

01 

8.200000 

OE 

01 

8.4CC000 

OE 

01 

8.60C000 

OE 

01 

8.80000C 

OE 

Cl 

9.000000 

OE 

01 

9. 200000 

OE 

01 

9.4CC000 

OE 

01 

9.6CC00C 

OE 

01 

9.800000 

OE 

Cl 

NUMBER  1 
MINI  MU  M=  2 


. 91 2092  6E 


01 


Lateral  pressure  p(y)  (kg./m, 


DfcPENDF  NT  VlA 

! 


I 

RI ABLE  (S 
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PLOTTING  SYMBOL 
DEPENDENT  V RI £BL  E 


Depth 

!ft 

.) 

I NDEPEND 

ENT 

VARIABL 

E 

0.0 

I. 000000 

OE 

00 

2.oooooq 

CE 

CO 

3.000000 

OE 

00 

4. cocood 

OE 

00 

5.0CC000 

OE 

00 

6.000000 

OE 

CO 

7.000000 

OE 

00 

8.000000 

OE 

00 

9.000000 

OE 

00 

1. cccood 

OE 

01 

l.  looood 

OE 

01 

1.200000 

OE 

01 

1.300000 

OE 

01 

1.40000C 

OE 

01 

1.5C00CC 

OE 

01 

1.60000C 

OE 

01 

1.700000 

OE 

01 

1.800000 

OE 

01 

1.900000 

OE 

01 

2.0CC00C 

OE 

01 

2. 10000C 

CE 

Cl 

2.200000 

OE 

01 

2.3C0000 

OE 

01 

2.400000 

OE 

01 

2. 50000C 

OE 

01 

2.60C00q 

OE 

01 

2.700000 

OE 

01 

2 • 800000j0E 

01 

2.9000000E 

01 

3.ooooocjoE 

3.  IOOOOOCE 

01 

01 

3.200000 

OE 

01 

3.300000 

OE 

01 

3.400000 

OE 

01 

3. 5C0000 

OE 

01 

3.600000 

CE 

01 

3.700000 

OE 

01 

3.800000 

OE 

01 

3.900000 

OE 

01 

4. CGCOCC 

OE 

01 

4.  10000C 

OE 

Cl 

4.20000C 

OE 

01 

4.300000 

OE 

01 

4.4COOOC 

OE 

01 

4.500000 

OE 

01 

4.60000C 

OE 

01 

4.700000 

OE 

01 

4.800000 

OE 

01 

4. 9C000C 

OE 

01 

5.000000 

OE 

01 

NUMB  ER 
MI  NI  MUM  - 

I 

* 

♦ 
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Depth  (ft.) 


Lateral  pressure  pCy) 


(kg. /m. ^ ) 


5 • 10000CCE  01 
5 • 2C00000E  01 

P S 3 Q99PQJQ?  2915678901  234567  8901  2345678901  2345678901  2345678901  23 


5.40C00CCE  01 
5.50000C0E  01 

5 • 6000000E  Cl  * 


5 ^70000© 

0E 

20156789 

0 123456789 

01  23456789 

01  23456789 

01  23456789 

01  2 3 4$  6 7 8 9 

5.800000 

OE 

01 

* 

5.9CC000 

OE 

01 

* 

6. CC000C 

OE 

01 

* 

6. 10000C 

CE 

Cl 

* 

6.200000 

OE 

01 

$ 

6.300000 

OE 

01 

* 

6.400000 

OE 

01 

* 

6. 500000 

OE 

01 

* 

6.600000 

OE 

Cl 

* 

6.700000 

OE 

01 

* 

6.800000 

OE 

01 

* 

6.900000 

OE 

01 

* 

7.0CC000 

OE 

01 

$ 

7.10000C 

OE 

01 

* 

7.200000 

OE 

01 

* 

7.300000 

OE 

01 

* 

7.4CC00C 

OE 

01 

* 

7.5CC00C 

OE 

01 

* 

7.600000 

OE 

Cl 

* 

7.700000 

OE 

01 

£ 

7.80000C 

OE 

01 

* 

7. 90COCC 

OE 

01 

* 

8. 0CC000 

OE 

01 

* 

8.100000 

OE 

Cl 

* 

8.200000 

OE 

01 

* 

8.300000 

OE 

01 

* 

8.400000 

OE 

01 

* 

8.50C00C 

OE 

01 

* 

8.600000 

OE 

Cl 

* 

8.700000 

OE 

01 

* 

8.800000 

OE 

01 

* 

8.9CC000 

OE 

01 

* 

9 • 0 0 CO  0 0 

OE 

01 

* 

9. 10000C 

CE 

Cl 

* 

9.200000 

OE 

01 

* 

9.300000 

OE 

01 

* 

9.4C000C 

OE 

01 

* 

9. 5COOCO 

OE 

01 

* 

9.600000 

CE 

01 

* 

9.700000 

OE 

01 

* 

9.800000 

OE 

01 

* 

9.9CC000 

OE 

01 

* 

TOTAL  NO. 

OF  INCREM 

ENTS  = 100 
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•#** 


(2)  Example  No.  2.  The  characteristic  functions 


are: 


Density- function: 

y = 45.0  lb/cu.ft. 

Ratio- function: 

2k 


k (y) 


o 


1+0 . 06y 


(constant) . 


where, 

k = 0.333. 
o 

Friction- function: 

X (y)  = y | cos  (y)  | 

where , 

y = 0.425. 


The  system  solved  is: 


aq  _ 

dy 


= y - H . ! cos (y)  I - 1+0,06y  * * 


2k 


q ( 0 ) = 0 
Geometry: 


H = 100.00  ft. 

R = 4.99  ft. 

Ay  = 1.00  ft. 

The  results  show  that  the  lateral  pressure  is 
larger  than  in  the  previous  case,  basically  because  of  the 
periodic  decrease  in  the  value  of  the  friction- function . 
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NUMERICAL  SOLUTION  OF  THE 
ANALYTIC  MODEL 

Example  No.  2 


DIMENSION  A ( 500 ) ,XX(500) 

EXTERNAL  FUN 

1 FGRMATOriO.  0,215) 

2 FORMAT  ( 1 P:  1 > 7 X , 4 4 H $ C L U T I 0 M OF  0 Y / DX = F UN ( X , Y ) BY  R.  K2  SUBROUTINE// 
1 1H, luX , 2HH= , F7d  3 ,2X,3HX0  = , F70  3t 2X ,3HY0= , F7o  3// 

" 2 1H , 12 X , IF  X ,1 8 X » 4H Y ( X ) / / ) 

3 FORM A T ( 1H , 1 0 X , F 5 . 2 , 10X , E 1 5. 8 ) 

“10  KEaD(5,1)  XDfYOt  H*  jNT»"lfcNT 

20  WR I T L ( 6 , 2 ) H,XJ,Y0 

CALL  RK2(FUN,H,  XU  , YO  , J N T , I E NT  , A ) ' 

STEP=FLGAT( JNT)*H 
X = X 0 

D0_  30  I = 1 , I E NT 
XX C I )=  X 

X=X+STEP  __  

A (I  ) = A ( I ) F ( (0.66  6 ) / ( 1 o 0 ♦ 0 "•  0 6 ^ XX  ( I ) ) ) 

B = A(  I )f40  P82  _ 

WRITE (5,751 ) X ,8 
751  FORM  AT ( 1 OX , FI Oo 5 , 1 OX , F 1 5 • 6 ) 

30  WRITE (6, 3)  X , A(  I ) 

CALL  A M A x ( A t I ENT, VALMAX , I SUB ) 

CALL  AMIN(A, I ENT, VALMIN, ISU3 ) 

CALL  PLOPd  ,IENT,1,VALMAX,VALMIN,XX,A)_ 

" CALL  PL 0 P ( U , I ENT , 1 , V A L M A X , V A L M I N , X X , A ) 

_ V A L M A X = 2 « 0 X • V A L M A X 

C A LL  PLOP (1,1 EN T , 1 , VALMAX ,VALMI N, XX, A) 

CALL  PLOP ( 0 , PENT , 1 , VALMAX , VAL MI N, XX ,A) 

WRITE (6,750)  I ENT 

7 50__F 0 RM A T ( * TOTAL  NO.  OF  INCREMENTS^*  ,15)  

40  contTnue 

__  _ CALL  EXIT  __ 

END 


^FUNCTION  FUN  ( X , Y ) 

FLA MO A = 0,425  * APS  ( COS (X)  ) 

_F  UN  = 4 5 o U — ( r L a M D A / 4 . 9 9 ) «•'  ( ( J » 6 6 6 ) / ( 1 « + ( 0 0 0 6 Y X ) ) ) ••  Y 
RETURN 
END  _ 
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SOLUTION  OF  D Y / 0 X = F IJ  N ( X,YI  BY  RK2  SUBROUTINE 


H=  loOOO  XO=  0.0  Y0=  0.0 


Depth  (ft.) 

Lateral  pressure  p(y) 
(kg./m.2  and  lb./sq.  ft.) 

Depth  (ft.) 

Lateral  pressure  p(y) 
(kg./m.2  and  lb.7sq.ft.) 

lo 00000 

143.296326 

"1  267  000  0 0 

1203.4553 22 

1.00 

0.29351 974 E 02 

26.00 

G.24753295E  03 

2.00000 

271. 109863 

1 27.00000 

1229.616455 

2.00 

0. 55532578E  02 

27.00 

0.251 86 740E  03 

3. 00000 

374.987793 

28.00000 

1 239,35^980 

3 o 00 

0.  76  8 1 03 03 F 02 

28.00 

0.25386218E  03 

4. 00000 

462.814453 

29. 00000 

1 2 4 2 o 2 0 507° 

4.00 

0. 94 3002 01 E 02 

29.00 

0.25444595^  03 

5. 00000 

551.970947 

30. 00000 

1253.906253 

o 

o 

. 

in 

0. 1 1 306248E  03 

30.00 

0.25786694E  0? 

60  00000 

622.187500 

31.00000 

1269.505659 

6.00 

0. 12744522E  03 

31.00 

0, 26003 809F  03 

7.00000 

673.049561 

32.00000 

1270.45922 3 

7.00 

0.1 388877 0E  03 

32.00 

0.26023340E  o3 

8.00000 

742.216064 

33.00000 

1283. 483  15-+ 

8.00 

0. 1 5203 1 1 6E  03 

33.00 

0.  2 62  901 12E  0 3 

9. 00000 

794.481689 

34. 00000 

1 295, 376953 

o 

o 

. 

C' 

0.16273497E  03 

34.  00 

0. 2653374 JE  03 

10. 00000 

830.998291 

35. 00000 

12  95.2  45  c 5 0 

1 0 . 00 

0. 17021 680E  03 

35.  00 

0.2653105 5 E 03 

11.00000 

878.069336 

36. 00000 

1 304  o 248  29  1 

11.00 

0. 17985855E  03 

36.00 

0.26715454E  03 

12.00000 

920. 593262 

3 7®  OOUOO 

1316. 739502 

12.00 

0. 18856891E  03 

37.00 

0.26971 31 3E  03 

13. 00000 

945.094971 

38.00000 

1316.284131 

13.00 

0. 19358768E  03 

38.00 

0.26960962E  03 

14. 00000 

973.622314 

39. G00U0 

1321.830076 

14.00 

0.  2004  552 3 E 03 

39.00 

0.27u7558oE  o3 

15. 00000 

1014.670654 

40, 00000 

1335.061279 

*— * 

. 

o 

o 

Q. 20783916E  03 

40.00 

0.2734660 6 E 03 

16. 00000 

1031.829346 

41.00000 

133^.742920 

16.  00 

0.  21  13  5385E  03 

41.00 

0. 27340J 8 BE  o3 

1 7. 00000 

1055.226318' 

42o  0000O 

1 33 7 o 577393 

17.00 

0.2161 4633E  03 

42,00 

0. 27398 145E  03 

18.00000 

1087.^48486 

43. 00000 

135io69l650 

18.00 

0o  22  2 74655  E 03 

43 , 00 

0,2 76372 56 E 03 

19. 00000 

1100. 272461 

44. OoOOO 

1351,996626 

19.00 

0. 22 53 7331 E 03 

44.00 

0 , 2 7 o 9 3 5 0 6 E o 3 

20. 00000 

1116.038574 

45. 00000 

1352, 6403^1 

2 0 o 0 0 

0. 22860275E  03 

45.00 

0 , 2 7 7 U 6 6 6 9 E 03 

21.00000 

1146.063477 

46. 00000 

1366,250739 

21.00 

0o  234752  90E  03 

46.00 

0,279854 74 E 03 

22. 000 00 

115b. 542725 

47. 00000 

1367, 549305 

22.00 

0. 23689941 E 03 

47.00 

0.28012085H  03 

23.00000 

1106.546143 

48.00000 

1 36  6.  55102  5 

23.00 

0. 23 8 94 846 F 03 

48.00 

0. 27991 626F  03 

24. 00000 

1193.109619 

40. 000 JO 

1377.490234 

24.00 

0. 2443 8 95 3 E 03 

49.00 

0,2821 56 98 E 03 

25. 00000 

1202. 682861 

5 o. ouu  30 

1330. 113525 

2 6. 00 

0o  246350 4 5E  03 

50.  00 

0.2826 74 34 F o3 
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Depth  (ft.) 

Lateral  pressure  p(y) 
(kg./m.2  and  lb./sq.  ft.) 

Depth  (ft.) 

Lateral  pressure  p(y) 
(kg./m.  and  lb./sq.  ft. 

51, 00000 

1378.020752 

76. 00 000 

1437.665527 

51.  O') 

0.28226562c  03 

76  o 00 

0.  294482  91 E 03 

5 2. 000 JO 

1 386. 647705 

77. 00000 

1442.021729 

52,00 

0*  284032  71 E 03 

77.0  0 

0.29537524E  03 

53, 00000 

1390. 816895 

78.00000 

1 445 o t>5  i 1 2 3 

53,00 

0. 2 648 86 72 F 03 

78,00 

0. 2 9611 865 E 03 

54,00000 

1388. 107666 

79. 00000 

1442.462891 

5^.00 

0.28433179c  03 

79  0 00 

0.29546558E  03 

55. 00000 

1394.753662 

80. 00000 

1445. 115967 

55.00 

u .285693128  03 

80.00 

0,296 0J903E  03 

56. 00000 

140U. 606934 

8 l.OCOOO 

1449.770264 

56.00 

(Jo  28o3920qE  0 3 

81.00 

Co  296962 40 E 03 

57.00000 

1397.689209 

82.00000 

1446. 801 270 

5 7.00 

0. 2 862 9443 E 03 

82.0  0 

Go  2963 542 5E  u3 

56.00000 

1401.959961 

83. ubOOu 

1 44  7 o 860840 

58,00 

0.28716919E  03 

83,00 

0.29657129E  03 

50. 00000 

1408. 879883 

94. 00 COO 

1453.459229 

59.00 

0.2885  866  7E  03 

8 4 o 0 0 

0.29771 6 02 E u3 

60. OOOUO 

1406. 1 39893 

85. 00 COO 

1450. 9 oO 93 7 

60.00 

0. 28802539E  03 

85.00 

0.2972U630E  03 

61. O'OOjO 

1408.285156 

86o  00* j'JO 

1 4 5 0 . o ? 0 4 3 • ) 

61.00 

0.28 846484E  03 

86.00 

0. 2971 5088 E 03 

62.00000 

1416. 230469 

87.00000 

1457.26733^ 

62.00 

0.29009229E  03 

8 7.00 

0. 298498 u5E  03 

63.00000 

1413. 994385 

83.00000 

1455.469707 

63.00 

0.28 96342 8E  03 

88.00 

0. 2981 2573 E 03 

64. 00  00 u 

1414.384277 

8 9.  00  000 

1454. 107666 

64.00 

C. 28971411E  03 

89.00 

0. 2978508 3E  03 

6 5.  oo  no 6 

1423.400879 

90, 00000 

1461.056396 

65.00 

0*291561 04E  u3 

90.00 

0.29927417E  03 

66. COOOO 

1421,937256 

9 1. OOOUO 

1460. 10 7 6 6 o 

66.00 

0.29126123E  03 

91.00 

0, 29907933 E 13 

67. 00000 

1420.922852 

92.00000 

1457. 953857 

67, 00 

0.29105347E  03 

92.00 

0. 29863 867E  03 

68.00000 

14  30©  07421 9 

93. OOOUO 

i 46  3.521240 

6 8.00 

0.  292  92  7 98  F 03 

93.00 

0. 299779 J5F  03 

69. OOUOO 

1429. 617676 

94. 00000 

1463.615479 

69.00 

0.29283447E  03 

94  „ 00 

0. 29979834 E 03 

70. 00000 

1427.547363 

95. 00000 

1 4 6 0 o 92  77  34 

70.0  0 

0.29241O40E  03 

95.00 

0. 299267 80E  03 

71.00000 

1434.890625 

° 6 o 00000 

1 '+65. 259033 

71.00 

0. 2939145 5E  03 

96.  0 0 

0.3001 3 50 IE  03 

72.00000 

1435. 658203 

97. 00000 

166  6.  5 2 7 1 u J 

72,00 

0.294071 78 E 03 

97.00 

0.  30  03  96-7 8E  0 3 

73.00000 

1432.895503 

9 3. oo  odd 

] 463. 6 o 42 C 9 

73.00 

0. 293505 °6E  03 

Q3 . 00 

0. 2 c978784  E 03 

74. 00000 

1 4 3 8 o 64160  2 

99.00000 

1466.6lj7^1 

74.00 

0 o 2 9 4 6 8 2 8 6 F 03 

99,00 

0.3GG45288E  o3 

7 5. 0 0 0 0 0 

1 440.  78930  7 

109, OOUOO 

14  69.  3^2529 

75,00 

0.29512230E  03 

M: 

0.30097 144 f 03 
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PkQTTING  symbol  _ __ . 

DEPEND E N T VARIABLE  NUMBER  1 

Depth  ( ft . ) Lateral  pressure 

‘"INDEPENDENT  M I N I MUM  = 2 . 93  5 1 9 7 A-  fc  01 

VAR  I ABL  E I J I | 


0.0 

2 • 0 0 0 0 0 0 C E 00 
4. 00 000 OLE  00 
6 » 000 000 L1  F 00 
8. 0000000 F 00 


* 


* 


* 


p(y)  (kg./m.2) 

DEPENDENT  VARIABLE ( 


o ) 

I, 


1. 00 000 00 E 01 
1 o 20 000 00 E 01 
1.4000000F  01 
,1. 6 0 000 GO  E 01 
1.800000OF  01 
2. 00 000  DOE  01 
2. 20 000 00 E 01 
2 • 40 00 000 E 01 


2. 6 0 00000 E 01  _ Y 

2. 800"6000E  01  ' ' ' ' v 


3. 00 000 00 E 01 
3 •2000 00 OE  01 
3. 4000  CO OE  0 1 
3.6000  0 00  E 01 

3.  8O0C000E  01 

4 .0  0 0 0 0 6 0 E 6"  T 
4. 20C000OE  01 
4. 400000° E 01 
4 . 6 0 C 0 0 0 0 E 0 1 

4.  8000000  p 01 
5.0000  0 0 0 E 01 
5. 200U000E  01 
5.40000 00E  01 
5. 6 0000 OOF  01 
5. 80 00000 r 01 
6. 00 000 00 E 01 
6. 200000OE_0i 
6 • 4 00  0 00  OF  ~ 0 1 


* 

* 


x- 

x- 

£ 

* 

3>i 


X 

A 


X 


A 


V 


6 

o 6 0 0 0 0 0 0 E 

0 

1 

6 

. 8000000 F 

0 

1 

7 

. 00 000 00 E 

0 

1 

7 

. 2 0 66 000 E 

0 

f 

7 

.40 000 00 E 

JJ 

1 

7 

o 6 0 6 G 0 00  E 

~0 

1 

7 

. 8 0 00  o 0 0 F 

0 

1 

B 

o 

O 

o 

o 

c 

r, 

O 

m 

0 

1 

8 

» 2 COCO OOF 

(1 

1 

8 

o 4 0 0 0 0 0 0 F 

0 

! 

8 

,.600pqpoF 

n 

1 

8 

• 80 00000 F 

0 

1 

9 

• 000C00QE 

0 

1 

9 

. 2000000E 

0 

1 

9 

. 400000 OE 

0 

1 

9 

.60  00  0 00  E 

0 

1 

9 

.800000  OF 

o 

1 

A 

V 


A 

A 

A 

x 

■*. 

A 


x 


A 


* 
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PLOTTING  SYMBOL  

DEPEND  EM T ' 'V AP I A B L E~ N UM3E R 1 

Depth  (ft.)  Lateral  pressure 

INDEPENDENT  ’ MINI  MIJM  = 2.  935 1 D74E  0 1 
VARIABLE  I I I I 

0.0  ' ‘ ‘ * 

1 « 0 0 0 0 0 0 0 E 00  * 

"2".  U 0 0 0 0 0 OE  0 0 ‘ ~ * ’ 

3.0000000F  00’  * 

A.  0000000 F.  00  ‘ * ” 

5._00  00000E_00  _ _ _ * _ _ _ 

6*  OOOOOOOE  00  " ‘ ~ * 


7 . 0 000 0 00 E_p  0 
8. OOOOOOOE  TO 
9.0000000F  00 
I . OOOOOOOE  7)1 
1. 1 0000 00 E 01 

I. 2 0000 00 E 01 

J. ±A9J>o  £? £ JIL 

i.4Gcooo  o e o f 
1 • 5C00000 E 01 
1.60U0000E  01 


p(y ) (kg. 

DEPENDE 


W700000JE  01 
1 .8 000000 E 01 
U 9000000 F 01 
2'.  OOOOOOOE'' OT 
2.1 0000  OOF  01 
2. 20000 OOE  01 


* 

- -■.4 


2.3000000E  01  ■* 

2. 4 0 CO 000 E 0 1 ' ' ‘ ' * 

2. 50 CC 00 OF  01  _ _ __  * 

2. 60000 ODE  oV  ""  ” ’ ‘ 

2. 7 000000 E 01  * 

2. 8 000 00 OF  01  * 

2. 9000000E  01  ■* 

3.  OOOOOOOE'  01  " " " * 


3 • 1 0 Of  ’ 00  0 F 01 
“3.  2 0000 00 E'  01 
3.30C0000E  01 
3.400  0 CO  0 E 0 1 
_3.  500 000 OE  01 
3.6 0000 OOE  01 
3.7000000E  01 

3.86oooo6f“'oT 

3 o 9GOOUOOF  01 
4. OOOOOOOE  oi 
4. 1 0000 OOE  01 
4.  2 0 000 OOE-  01 
4*_3000090F  01 
4. 40 JO 000 E 01 
4, 5 000 000 E ul 
4.60  LOO OOF  ul 
4 • 7 U 0 0 0 0 0 E 01 
4. 8 0 00 000 E 01 
4.9O00000F  01 
”5.  OCOC'OOOE  01 


■a 

v 

* 

* 
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Depth  (ft») 


Lateral  pressure  p(y)  (kg./m.2) 


5.1.O0000UF  til 

5.  ?doooooi?  oi 


5*  30000 00E  01 
5 • 400000 JE  01 
5 • 5 OOQOoOE  01 
5#  6 0 000 00 E 01 
5. 70 000 00 E 01 
5 • 8 0 0 0 0 0 0 E 01 
5 • 9 0 00 0 J 0 E 0 1 
6. 0000000 E 01 
6« 1000000 F 01 
2 OuuOOOE  01 
6. 3 000 00 OF  01 

6.  400000 J E 01 

6o5 C 00 UOO E 4 1 

6. 60 000 00 E 01 
6_?7q0v''0u0F  01 

6*.  a 00 o’ 000 E 01 
6. 9000n00E  01 
7.00000  .10  E "oi 

7. JOCOOOOE  01 
7 • 2 0 0 0 0 0 0 F 01 


7.  3 0^00 0 00 E_  01 
7*4000000E  01 
7c5000000F_01 
7. 6 000 dob E 01 
7.70000UOE  01 
7.8bobbbo~E  01 
7 •90  0000  Uj  01 

8.  “bob 00 OOF  01 
8* 10 00 000 E 01 


8. 

2 0 000 

00  E 

91 

8 o 

3001)0 

OOF 

01 

8 c 

4 0 000 

90  E 

01 

8. 

5 0 0 ( 1 0 

OOF 

01 

~8o 

6 0 COO 

OOE 

01 

8 o 

70000 

00  E 

01 

8. 

8000  0 

00  E 

01 

8. 

90000 

OOE 

01 

9. 

ob  obb 

0 JE 

0 1 

90 

10  000 

90  c 

01 

90 

2 0000 

OOE 

01 

9. 

3 0009 

0 0 F 

01 

9. 

4 000  0 

OOE 

01 

9# 

5 0 ( j 0 9 

90  E 

01 

9 o 

60009 

00  E 

01 

9o 

70  000 

OOF 

91 

9 o 

8 000  0 

OOE 

01 

9. 

9000U 

OOE 

01 

TOTAL  NCI.  OF  INC  RE  ME  NT  S = 


¥ 

V 
* 
¥ 
¥ 
■¥ 
¥ 
* 

Y 


¥ 

Y 

Y 
¥ 
Y- 

Y 
¥ 


¥ 

Y 

Y 

Y 
¥ 

Y 
¥ 


Y 
* 
A 

Y 

Y 
¥ 
4c 
¥ 

¥ 

¥ 

¥ 

¥ 


¥ 


•V 

¥ 

¥ 

5 


¥ 

¥ 

•* 


10  0 
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(3)  Example  No.  3.  This  example  is  identical 


with  example  No.  2,  except  for  the  value  of  the  friction 
coefficient  y,  i.e.,  y = 0.850.  The  large  value  assumed 
for  y and  the  periodic  behavior  of  \ (y)  produce  pressure 
fluctuations  in  the  solution.  For  details,  see  the  plots 
accompanying  this  example. 
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NUMERICAL  SOLUTION  OF  THE 
ANALYTIC  MODEL 

Example  No.  3 


DIMENSION  A ( 500 ) , XX ( 500 ) 

EXTERNAL  FUN 

1 F0RMAT(3F10.0,2I5) 

2 FORMAT  (1H1,7X,44HSOLUTION  OF  DY/DX=FUN  ( X , V ) BY  RK2  SUBROUTINE//* 
11H, 10X,2HH=,F7. 3 , 2 X , 3HX0= , F7. 3 , 2 X , 3HY0=  , F7. 3 / / 

21H, 12  X, 1HX,18X,4HY(X) //) 

3 FORMAT (1H,10X,F5.2,10X,EI5.8) 

10  RE AD ( 5 , 1 ) XO,YO,H, JNT, IENT 

20  WRITE(6,2)  H,XO,YO 

CALL  RK2(FUN,H, XO, YO, JNT, IENT, A) 

STEP=FLOAT(JNT)*H 

x=xo 

DO  30  1=1, IENT 
XX  C I )=X 
X=X  + STEP 

A(I)=A(I)*( (0.666) /(  1 .0  + 0. 06*XX ( I ) ) ) 

B=A ( I )*4.882 
WRITE(6,751)  X , B 
751  FORMAT ( 10X,F10. 5 ,10X,F15.6) 

30  W R I T E ( 6 , 3 ) X , A ( I ) 

CALL  AMAX(A, I ENT ,VALMAX,  ISUB) 

CALL  AMIN( A, I ENT, VALMIN, ISUB ) 

CALL  PL  OP ( 1 , 1 ENT,1 , VALMAX, VALMIN, XX, A) 

CALL  PLOP (0, I ENT,1, VALMAX, VALMIN, XX, A) 

VALMAX=2.0*VALMAX 

CALL  P L 0 P ( 1 , IENT, 1 , VALMAX, VALMIN, XX, A ) 

CALL  PLOP (0, I ENT ,1, VALMAX, VALMIN, XX, A) 

WRITE (6,750)  IENT 

750  FORMAT ( • TOTAL  NO.  OF  I NCRE MENT S = • , I 5 ) 

40  CONTINUE 
CALL  EXIT 
END 


FUNCTION  F UN ( X , Y ) 

F L AMD A = 0.850  * ABS  ( COS(X)  ) 

FUN  = 45.0  - (FLAMDA/4.99)*( (0.666) /(l.+(0«06*X) ) )*Y 

RETURN 

END 
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SOLUTION  OF  DY/DX=FUN( X, Y)  BY  RK2  SUBROUTINE 


H=  1*000  X0=  0.0  Y 0=  0.0 


Depth  ( f t . ) 

Lateral  pressure  p(y) 
(kg./m.2  and  lb./sq.  ft.) 

Depth  (ft.)  (kg 

Lateral  pressure  p(y) 
./m.2  and  lb./sq.  ft.) 

1.00000 

' 140.367920 

26. 00000 

982.937500 

1.00 

0.28752136E  02 

26.00 

0. 20 1339 1 1 E 03 

2.00000 

266.310791 

27. 00000 

1004.559814 

2.00 

0.  5454953  OE  02 

27.00 

0. 2057681 1 E 03 

3.00000 

358.965576 

28.00000 

10U8. 114502 

3.00 

0. 7352842 7E  02 

'28.00 

0. 2064962 5E  03 

4. 00000 

432.407959 

29.00000 

1001.284912 

4.00 

0. 88571915E  02 

29.00 

0. 20509731 E 03 

5.00000 

517.428223 

30. 00000 

1017.449219 

5.00 

0. 10598697E  03 

30.00 

0.20840833E  03 

6. 00000 

574.696533 

31. 00000 

1024.171631 

6.00 

0. 11771 745  E 03 

31.00 

0.20978531E  03 

7.00000 

612.612305 

32.00000 

1015.980713 

7.00 

0. 12548390E  03 

32.00 

0.20810753E  03 

8.00000 

671.320557 

33.00000 

1027.667725 

8.00 

0. 13750935E  03 

33.00 

0. 2 1050142F  03 

9. 00000 

713. 176514 

34.00000 

1037. 849365 

9.00 

0.  146082  90E  03 

34.00 

0.21258694E  03 

10. 00000 

731.459717 

35.00000 

1029.277100 

10.00 

0.14982790E  03 

35.00 

0. 2 1083107E  03 

1 1.00000 

772.277344 

36. 00000 

1035.721924 

11.00 

0. 158188 75 E 03 

36.00 

0. 21215120E  03 

12.00000 

808.063721 

37.00000 

1047.955566 

12.00 

0.16551903E  03 

37.00 

0. 21465706E  03 

13.00000 

815.428223 

38. 00000 

1039. 910156 

13.00 

0. 16702754E  03 

38.00 

0.21 300908E  03 

14.00000 

841.225830 

39.00000 

1041.942871 

14.00 

0. 17231 1 77E  03 

39.00 

0.21342546E  03 

15. 00000 

873.391113 

40.00000 

1056.232178 

15.00 

0.17890031E  03 

40.00 

0. 21635239E  03 

16.00000 

875.291748 

41.00000 

1049.401123 

16.00 

0.17928960E  03 

41.00 

0.21495316E  03 

17.00000 

889.916992 

42.00000 

1047.908936 

17.00 

0. 18228 53 7E  03 

42.00 

0. 2 1464749E  03 

18. 00000 

920.904297 

43. 00000 

1064.295410 

18.00 

0.18863260E  03 

43.00 

0.21 800401 E 03 

19. 00000 

920.644043 

44. 00000 

1059. 1 84326 

19.00 

0. 18857933E  03 

44.00 

0.  21695  709E  03 

20.00000 

926.958984 

45. 00000 

1054. 938477 

20.00 

0. 18987282E  03 

45.00 

0. 2 1608737E  03 

21. 00000 

958.094727 

46. 00000 

1071.118403 

21.00 

0. 1 9625046E  03 

46.00 

0.21940157E  03 

22.00000 

957.71826 2 

47.00000 

1068. 181 152 

22.00 

0.19617337E  03 

47.00 

0.21879996E  03 

2 3.00000 

957.856445 

48.00000 

1061.945063 

23.00 

0.19620168E  03 

48.00 

0.21 7522  55E  03 

24.  000  JO 

986.1 10107 

49.00000 

1074.458252 

24,00 

0.20198897E  03 

49.00 

0.  22008569E  03 

25. 00000 

987.125000 

50. 00000 

1074. 126953 

25.00 

0.202196 85 E 03 

50.00 

0. 22001 785E  03 
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Depth  (ft.)  (kg 

51. 00000 

51.00 

52. 00000 

52.00 

53.00000 

53.00 

54. 00000 

54.0  0 

55. 00000 

55.00 

56. 00000 

56.00 

57.00000 

57.00 

58.00000 

58.00 

59.00000 

59.00 

60.00000 

60.00 
61.00000 

61.00 
62.  d)0000 
62.00 
63.00000 

63.00 

64.00000 
64.  0(J) 

65.00000 

65.00 

66.00000 

66.00 

67. 00000 

67.00 

‘ 68.00000 

68.00 

69.00000 

69.00 

70.00000 

70.00 

71.00000 

71.00 

72.00000 

72.00 

73.00000 

73.00 

74. 00000 

74.00 

75. 00000 

75.00 


Lateral  pressure  p(y) 
./m.  and  lb./sq.  ft.) 


Depth  (ft,) 


1066.703369 
0.21849722E  03 

1076.064209 
0. 22041 464E  03 

1078.594971 
0. 22093306E  03 

1070.662598 
0. 21930823E  03 

1077.346924 
0. 22067741 E 03 

1082.862061 
0.221 80 708E  03 

. 1074.988281 

0. 22019427E  03 

1078. 346436 
0. 220882 13E  03 

1085.831299 
0. 22241 528E  03 

1078. 565918 
0. 2209271 1 E 03 

1078. 958008 
0.22100 739E  03 

1088.293457 
0. 22291 963E  03 

1082.092285 
0. 2 2 1 64941 E 03 

1080.061035 
0. 22 1 23334E  03 

1091.265381 
0.22352837E  03 

1086.489014 
0.22254999E  03 

1082.536865 
0. 22 1 74049E  03 

1094.  144043" 
0.2241180GE  03 

1091.114746 
0. 22349751 E 03 

1085.745850 
0.22239777E  03 

1094. 772705 
0.22424676E  03 

1093.807617 
0. 224049 1 0E  03 

1087.571045 
0.22277162E  03 

1094©  332275 
0.22415657E  03 

1095. 642578 
0. 22442496E  03 


76.00000 

76.00 

77.00000 

77.00 

78.00000 

78.00 

79. 00000 

79.00 

80.00000 

80.00 
81.00000 

81.00 

82.00000 

82.00 

83.00000 

83.00 

84. 00000 

84.00 

85.00000 

85.00 

86. 00000 

86.00 

87. 00000 

87.00 

88. 00000 

88.00 

89.00000 

89.00 

90. 00000 

90.00 

91.00000 

91.00 

92.00000 

92.00 

93.00000 

93.00 

94. 00000 

94.00 

95. 00000 

95.00 

96. 00000 

96.00 

97.00000 

97.00 

98. 00000 

98.00 

99. 00000 

99.00 

100. 00000 


Lateral  pressure  p(y) 
./m,2  arid  lb./sq.  ft.) 

1089.028320 
0. 2230701 3E  03 

1093.831543 
0. 224053 99E  03 

1097. 539062 
0. 22481 343E  03 

1090.976807 
0.22346924E  03 

1093. 332764 
0.223951R6E  03 

1098.727295 
0.22505682E  03 

1092.633301 
0.22380858E  03 

1092.704590 
0.22382 3 18E  03 

1099.642090 
0. 22 524420E  03 

1094. 379639 
0.22416629E  03 

1092.557373 
0. 223 79300E  03 

1101.069092 
0. 225 5365 OE  03 

1096.937983 
0.22469034E  03 

1093.598145 
0.  2 240062 1 E 03 

1102.756104 
0.22588205E  03 

1100.022705 
0. 225322 17E  03 

1095.542725 
0. 22440454E  03 

1102.700684 
0.  22587071E  03 

1101.627197 
0. 22565083 E 03 

1096. 427979 
0. 22458585E  03 

1101.804932 
0.22568724E  03 

1102.580078 
0. 22584601 E 03 

1097.046875 
0.22471262E  03 

1100.867920 
0.2254953 2E  03 

1103.613037 
0.22605762E  03 
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PLOTTING  SYMBOL  _ _ * 

DEPENDENT  VARIABLE  NUMBER  1 

Depth  (ft.)  Lateral  pressure  p(y)  (kg./m.  ) 

“INDEPENDENT  MINIMUM*  2.8752136E  Ol"”  DEPENDENT  VARIABLE  Q. 

VARIABLE  I 1 I I II 

0.0  * 


I •jOOOOOOOE  00 

2# odooouoE  bo 

3.0000000E  00 
4«  OOOOOOOE  00 
5.0000000E  00 
6 • OOOOOOOE  00 
7.0000000E  00 
8* OOOOOOOE  00 
9. OOOOOOOE  00 
1#  OOOOOOOE  01 
1 • 1000000E  01 
1.2000000E  01 
1.3000000E  01 
1.4000000E  01 
1.5000000E  01 
1 • 6000000 E 01 
1.7000000E  01 
1 • 8000000E  01 
_1  • 9000000E  01 
2. OOOOOOOE  01 
2.1000000E  01 
2*  2000000E  01 
2.3000000E  01 
2. 40 00 000 E 01 
2.5000000E_01 
2.6000000E  01 
2.7000000E  01 
2.8000000E  01 
2*  9000000E  01 
3# OOOOOOOE  01 
3* 1000000E  01 
3«  20000  00E  01 
3.3000000E  01 
3.400CC00E  01 
3«  5000000E  01 
3«  6000000E  01 
3.7000000E  01 
3 • 8000000E01 
3.9000000E  01 
4. OOOOOOOE  01 
4. 1000000E  01 
4.2000000E  01 
4.3000000E  01 
4*  40 00 000 E 01 
4*  5000000E  01 
4«  6000000 E 01 
4* 7000000E  01 
4«  8000000E  01 
4«  90C0000 E 01 
5. OOOOOOOE  01 
5# 1 OOOOOOE  01 
5.2000000E  01 
5 • 3000000E  01 
5.4000000E  01 
5«  5000000E  01 
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Depth  ( f t . ) 


Lateral  pressure  p(y)  (kg./m*2) 


5.6000000E  01 
5.7000000E  01 
5.8000000E  01 
5.9000000E  01 
6*  0000000E  01 
6.1000000E  01 
6.  2000000E  01 
6*  3000000E  01 
6.4000000E  01 
6.5000000E  01 
6 • 6000000E  01 
6.7000000E  01  ___ 

6*  8000000E  01 
6.9000000E  01 
7.0000000E  01 
7* 1 000000E  01 
7.2000000E  01 
7.3000000E  01 
7*  4000000E  01 
7*  5000000 E 01 
7.6000000E  01 
7.7000000E  01 
7*  8000000E  01 

7 • 9000000E  01  _ 

8* 0000000E  01 

8. 1000000 E 01 

8*  2000000E  01 

8 • 300G000E  01 

8.4000000E  01 

8.5QGOOOOE  01 

' 8#  6000000 E 01 
8*  7000000 E 01 
8*  8000000E  01 
8.9000000E  01 
9*  0000000 E 01 
9*  1 000000 E 01 
9.2000000E  01 
9. 3000000E  01 
9.4000000E  01 
9*  5000000E  01 
9.6000000E  01 
9#  7000000E  01 
9.8000000E  01 
9.9000000E  01 

TOTAL  NO.  OF  I NC  RE  MENT  S=  100 


i 
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5.2 


COMPUTERIZED  GRAPHIC  REPRESENTATION  OF  THE 


ANALYTIC  MODEL 


The  models  represented  in  the  form  of  Eq.  (3.18) 
describe  three-dimensional  surfaces  over  the  y-q  plane. 
Graphic  visualization  of  these  surfaces  may  aid  in  the 
model  building  process.  Surfaces  describing  various 
analytic  models  were  plotted  automatically  by  a Calcomp' 
plotter.  The  program  listings  are  given  in  Appendix  B. 


5.2.1 


are: 


where , 


Examples  and  Results 

(1)  Example  No.  1.  The  characteristic  functions 


Density- function: 

Y(y)  = y0  + Sy 


Y = 45.0  lb/cu.ft. 
o 

<5  =0.06 

Ratio- function: 

2k 


k (y)  = 


yk 

l+(-^)y 


where , 

k = 0.333 
o 

y = 0.425 
R = 5.0  ft. 
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Friction- function: 


X = 0.425  (constant). 

The  surface  plotted  is: 

= Y(y)  “ | * k(y)q(y) 

0 £ y £ 120  ft. 

0 _<  q 400  lb/sq.ft. 

The  plotted  surface  can  be  given  in  different 
orientations,  registered  by  the  coordinate  system  at  the 
upper  right  corner  of  each  plot.  Three  representative 
plots  of  this  surface  are  given  in  Figs.  5.1,  5.2,  5.3. 

(2)  Example  No.  2.  The  model  is  the  same  as 
in  example  No.  1,  except  that 

X (y)  = y | cos  (y)  | . 

Three  orientations  of  the  resulting  surfaces  are 
given  in  Figs.  5.4,  5.5,  5.6. 

(3)  Example  No.  3.  The  model  is  the  same  as 
in  example  No.  1,  except  that 

X (y)  = y | sin  (y)  | . 

Three  orientations  of  the  resulting  surfaces  are 
given  in  Figs.  5.7,  5.8,  5.9. 
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RLPHR- 
BETR- 
G R M M R - 


-■o.  00000 
20.00000 
=45.00000 


XMIN=  0.00000 
YM  I N = 0.00000 
Z M I N - 22.35601 
slcck 


XMRX--  120.00000 
YMRX-  400.00000 
ZMRX-  52.  10335 


Figure  5.1 
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? 


XMRX  = 120.00000 
YMRX=  lj 00 . 00000 
ZMRX  = 02, 19000 


RL PH R -0.00000 
BETA=  MO. 00000 
S A M M R - LJ 0 . 00000 


XMIN=  0.00000 
YM  I N = 0.00000 
Z M I N - 22.30601 


Figure  5.2 
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ALPHA=0. 00000 
BETA-  160.00000 
G AMM A = 4 5 . 00000 


XM  I N = 0.00000 
YM  I N = 0.00000 
Z M I N=  22.30601 

BLOCK  5 


XM  AX  = 120.00000 
YM  AX  = 400.00000 
Z M AX  = 52.  19995 


Figure  5.3 
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Z - 45. 

+ .06  « X 

- . 056 

* ABS  (COS  (X)  ) 

« T / 

(1  + .0233  * 

ALPHA-0 . 

00000 

XMIN- 

0.00000 

XMAX- 

120.00000 

BETA-  20 

. 00000 

YMIM-- 

0. 00000 

YMAX- 

400.00000 

GAMMA- 4 5 

. 00000 

Z M I M = 

22.  .35601 

Zmax- 

52.  19995 

eiccK  i 

Figure  5.4 
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Z « 45 . + . Q6  * 
RLPHR=0 . 00000 
BE  T A = GO. 00000 
G R M M A --  4 5 . 00000 


. Q5fi  * RB5  ICO  j IX)  ) 
XM  I N = O.OCOCO 
YM  I N = 0.00000 
Z M I M=  22.35601 

ELCCK  2 


* Y / (1  + 0283  * X) 

X M A X - 120.00000 
YMAX-  400.00000 
ZMAX=  52 . 1 0930 


Figure  5.5 
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X 


z 


Z = 4 5.  + .06  * X 
ALPHA-0 . 00QQ0 
BPT  A-  180.00000 
GAMMA -40. 00000 


- . 006  * ABS  (COS  (X)  ) 
XMIN-  O.COOGO 
YM  I N = 0.00000 
Z M I N=  22.35601 

BLGCK  ri 


* Y / (1  + .0283  k X) 

XM AX=  120.00000 
YMAX-  400.00000 
Z M A X - 02 . 19990 


Figure  5.6 
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Z = 45.  + .06  x 
ALPHA-0. 0GG0Q 
BET  A = 20.00000 
GAMMA-45. QGOQQ 


. 056  v AB5  (SIN  (XI  I 
X M I N = 0.00000 
X M I N = 0.00000 
Z M I N = 25.42896 

BLCCK  1 


x T / (1  + .0283  X XI 

XMAX=  120.00000 
XM  AX  = 400.00000 
Z M A X - 52.  19995 


Figure  5.7 
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Z - 45.  + .06  V X - 
RlPHA-0 . 00000 
BETA-  100.00000 
GAMMA-45 . 00000 


. 056  x ABS  (SIN  (X) I 
X M I N - 0.00000 
Y M I N - 0.00000 
Z M I N = 25.42806 


x Y / (1  + .0283  x X) 

XM AX=  120.00000 
YM AX=  400.00000 
ZMAX-  52. 19995 


Figure  5.8 
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Z - 45.  + .06  * X 

- .056 

« ABS  (SIN  (X)  ) 

« T / 

(1  + .0283  x X) 

ALPH A=0 . 00000 

XM  IN  = 

0.00000 

X M fl  x -■ 

120.00000 

BFT  A = 140.00000 

TM  I N = 

0.00000 

YMflX  = 

400.00000 

GflMNA=4S. 00000 

ZM  I N = 

25.42896 

eLCCK  4 

ZMflX  = 

52.  19995 

Figure  5.S 
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(4)  Example  No.  4.  The  model  is  the  same  as  in 
example  No.  3,  except  that  the  ranges  were  extended  to  the 
negative  directions,  as  follows, 

”120  <■  y £ 120  ft. 

”400  _<  q £ 400  lb/sq.ft. 

Three  orientations  of  the  resulting  surfaces  are 
given  in  Figs.  5.10,  5.11,  5.12. 
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z 


Z = 45.  + .06  * X - 
RLPH A-0 . 00000 
BETA--  20.00000 
GAMMA --45  . QOQQQ 


. 056  * ABS  (SIN  (X)  ) 
XM  I N = -120. QOQQQ 
YM  I N=  -MOO . 00000 
Z M I N-  -335.28843 
ai ecK  1 


k Y / (1  + .0283  x X] 

XM AX=  120.00000 
YMAX=  400.00000 
ZMAX-  421.25616 


Figure  5.10 
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i = 45.  + .06  « X - 
ALPHA=0. 00000 
B6TA=  100.00000 
GA1MA=45 . 00000 


.056  x AES  (SIN  (X) I 
XHIN=  -120.00000 
YN  I N = -1400.00000 
ZMIN=  -335.28349 

F.LOCK  3 


* T / (1  + .0233  X X) 

XMRX  = 120.00000 
YMRX-  400.00000 
ZMflX-  421.25616 


Figure  5.11 
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Z = 45 
ALPHA- 
BET  A- 
GAMMA- 


Z 


* 


. + . QG  x 

0. QQQQQ 
180. QQQQQ 
45.00000 


. 056  x ABS  (SIN  (XI  3 
X M I N - -120. 00000 
YM  I N = -400.00000 
ZM  I N-  -335. 28848 

2LGCK  5 


x Y / (1  + .0283  x XI 

XM  AX-  120.0.0000 
YM AX - 400.00000 
Z M A X - 42  1.25616 


Figure  5.12 
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5.3 


DIGITAL  SIMULATIONS  BASED  ON  THE  ALGEBRAIC  MODEL 


The  computer  models  of  this  category  are  based  on 
the  algebraic  model  given  by  Eq.  (3.28).  This  equation 
is  programmed  in  subroutine  SPH.  Various  numerical 
experiments  were  carried  out  and  some  examples  and  results 
follow. 

5.3.1  Examples  and  Results 

(1)  Example  No.  1.  Janssen-type  (static)  conditions. 
The  characteristic  functions  are: 

Density- function: 

3 

y = 750  kg/m  (constant) . 

Ratio- function: 

k = 0.333  (constant) 

Friction- function: 

A = y = 0.450  (constant).  Normal  friction. 

2 

The  results  give  the  lateral  pressure  (kg/m  ) in 
tabular  form  for  99  m at  1.0  m intervals. 

(2)  Example  No,  2.  Same  as  example  No.  1,  but 
with  low  friction,  i.e.,  y = 0.225.  The  results  show  that 
the  lateral  pressure  increase  significantly. 
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DIGITAL  SIMULATION 
Example  No.  1 


Cl  MENS ION  ON S 
RE  A 3 11  1 1 Cl ) I 
lOY  FORMAT  (I  3 » F 5 . 
CALL  SDNSTYCD 
CALL  SRA  T I C{ 
CAL-  SWALLtWA 
CAL-  SPHIDNST 
CALL  OLTPTHP 
CALL  E XI  T 
END 


TY  ( 100  ) »RA[T  I C(  100)  ,W 
H f R HO 
2) 

NST Y, I H ) _ 

RATI C, IH  ) 

LL,IH)__ 

Y #RAT  10  tWA 
rtf  IH) 


ALL  (ICO)  , PH(IOO) 


LLflHt  PH  , RTC) 


__  _ 

SUB 

ROUTINE  SC 

n’styI cnsty 

tIH) 

DI  M 

i NS  I ON  DNS 

I Y C 1 0.0  JL 

DO 

25  J=  1 1 1 H 

. t 

25  ON  S 

tYC  J ) =7  50  . 

D 

RET 

URN 

END 

25 


SUB 

ELM 

CO 

HA 

RET 

END 


ROUTINE  SRlATICtRATIO 
ENS  ION  RAT|I  C(  IOC) 

25  J = 1 t IH 
01  J 1=0. 3 3 
JRN 


T 1 


t IHJ 




SUE 

DIM 

DO 

_ 25  HAL 
RET 
ENC 

ROUTINE  S V» 
b NS  ION  >'AL 
25  J=l, IH 
L ( J ) = 0 . A 50 
CRN 

ALL(  WA  LL  ,1 
- ( 100  ) 

rt  ) 
95 

- - — - - — 

SUBROUTINE:  SP 
diHens.icn  DNSfT 
cjt  n 1100) 


Hi DNSTY  ,RA 
Y (1  CU  , R A 


CO.  BLO  _K  “ 1 » I H 
5 C XMATRX  lK)=1.0f-RATIC(  K) 
.DO  95__IZ=2  , I F 


DET  ( 11  = 1-0 


IZP 

{ DO 

0 12  3 4 5 6 J-rK 


1=1  Z— 1_ 

|7  5 K ^ 1 f I Z W1 

(El  2 3 4 5 6 7 8 9 


0 1 23456789 

IMlhl-1 
IZPK^I  Z-K 
75  DETt  I)=DET  ( I M 1 ) *XMATRXE 
SUP=DNSTY{ IZ) 

CO.  Qb  K = 1 , I Z P 1 
Kl  = < +1 
IZM<=i  Z-K 
85  SUM=SJK4-CNSTYi 
FH(  I.Z)  - RA  TI  0 (|I 
95  CONTINUE 
RETURN 
ENC 


II C > W ALL  , 1 
TIO(IOQ.)jW 


r/A  LL  l K)  / FH 


IZMK  i *DET[E 
Z)*$UM 


i,  PH,  RHO) 
ALL  ( TOO  ) ,P 


3 


%) 

01  23456789  Loi  23456789 


IZMK  ) 


KI) 


Hi  100),  XMA 


fRX  E100 ) , 


0123456789 


0 1 2 3 4 5 6 7 8 


..SUB 

D 1M 
„ V»R 
200  FOR 
_ F£  T 
END 


ROUTINE.  OU 
ENS  ION  PH l 
ITEL2 ,200  ) 

^ AT  ( T 3 0 , F 1 
URN 


TPTl (PH, I H) 
100) 

1 P H ( P ) , P=2 
0.5  ) 


,1  H ) 


96 


Lateral  pressure  p(yl 
(kg. /m. ^ ) 


49  1 
_ 7 26 
954 
1175 
1390 


15 


96 


1800 

1995 

2135 

2369 

2548 

2721 


2 6 
30 
32 


6 5 
52 
10 


3353 

3512 

3656 

3796 

3532 


40 


64 


4191 


43 


4435 
4 5 52 
4665 


47 

46 

49 


50B4 

5181 


5.2 

53 


0 1 2 3 4 5 6 7 


5*55 


15 


75 

81 

84 


75 

£7 


.99576 

.97485 

.89380 

.96802 

.40332 

.3  581.9 

.14615 

.83£€7_ 

.65282 

.76660 

.35229 

.57471  _ 

.55521 

. 57CC7_ 

• 6 5 C 1 5 
.98  264 
.7 12 16 
.,97412 
. 90430 
. 63155 
.28365 

• 5 8 C4  7 
.84275 
.98828_ 
.52344 
.55859 
.20312 

. 55C78 

• 7 C 703 
.76562 
.61641 

• 95  703„ 

. 26562 

• £359*8  9 


0 1 2 3 4 5 6 


55 

56 

57 

57 

58 
55 
60 
60 
6 1 
62 
62 
63 

63 

64 

65 
65 


41.74605 

25.07422 

05  .90234 
84.30078 

6 C.  34766 
34. 1C537 
0 5.  65^  25 
76.05465 
42.27  109 

07  .660  16 
7C. 99609 
32.42578 
92. Cl  172 
49. 6C465 
05. 86715 
60 .24  2 19 


Lateral  pressure  p(y) 
(kg. /m. 2 ) 


66 

12.98828 

66 

64.14453 

67 

1 3 • 76953__ 

67 

61. 50234 

68 

06. 5 8584 

68 

52.875CC 

68 

97  .80078 

69 

40.41016 

69 

81.73437 

- 7 C 

21.62031 

7C 

6 C . 7 C 703 

70 

98.42187 

71 

35.00C0C 

71 

70.48437 

7,2 

04  .90234  _ 

72 

38.28906 

.L  _ 72 

7C. 66797 

73 

02.C7812 

73 

32. 54257 

72 

62.05375 

_ J 13 

90.7578  1 _ 

74 

18 .56  250 

_ 4 74 

4 5. 53 12  5_ 

*'•:  74 

71.60750 

74 

57.C585.5„ 

75 

21.67187 

75 

45.54297 

75 

58  .69  531 

15 

91  ,.15625 

- 76 

,12.94141 

U _ 76 

[3  4.  C 7.03.1 __ 

It 

54. 56641 

_ _ 16 

74.44  53  1 _ 

76 

93. 73C47 

17 

12.43259 

77 

30 .57422 

77 

48.17  187 

77 

65.24219 

77 

81.  75687.. 

" 7 7 

97.  85537 

_ „ 78 

12.43355- 

2 78 

28 .54297 

456778 

4 3»19922a 9 o 

_ 78 

57.41406 

78 

71.20312 

78 

84. 57812 

78 

97. 55076 

r*  79 

10 . 1328  1 
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DIGITAL  SIMULATION 
Example  No.  2 


LIM 
RE  A 
1C 1 FOR 
CA  L 
CAL 
CAL 
CAL 
CALL 
CAL 
END 


ENS  ION  DIMS 
D (1  t 1 01 ) I 
MAT  (I  3,F5 
L SD NST  V ( D 
SRATI  C(  R 
S WAL L(  WA 
S P H ( DN  S T 
0 LTPT1 (P 
E >1  T 


L 


TY(  100  ) ,RAtTI  C(  IOC)  ,W  ALL  (100)  ♦ 
H,RHO 
2) 

MSTY, IH) 

I G, 1H  ) 

L L * I H ) 

Y »RA  T 1 0 » WA[LL  * I H f PH,  R^O) 

H ,1  hi 


SUE 
DI  M 
DO 

25  .CNS 
RET 
• ENC 


ROUTINE  SD 
E NS  I ON  DNS 
25  J=lt  IH 
[TY(  J ) = 7 50  . 
URN) 


\|  STY<  DNSTY 
TY(IOO) 

0 


IH) 


SUB 
Cl.M 
CO 

25  PATH 
PE  r 
END 


< GU  T I NE  SR 
ENS  ION  RATI 
25  J = L , IH 
0{ J J=0.33 
URN 


AT  I 0 ( R A 11  0 
0I.10C). 


, I H) 


P H ( 100  ) 
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SUE 

DIM 

ROUTINE  SVs 
b NS  I ON  Vs'  AL 

ALL! WALL, I 
- ( 100) 

H ) 

— — — — 

DC 

25  VsAL 

25  J=N  IH 
L(  J ) = 0 . 22  5 

RET 
EN  C 

URN 

SUBROUTINE  SP 
__  _ C IM|bNS..  IU  N D NS 
D^T(  110  0(j 
“ E»C.  K=1  , IH 
5 C >MATRX  (K)  = 1.0 

S 5 I 1-2 i I H 

ft  i) -l.o 

1=12-1 


1 
i 

0 1 2 3 4 5 


7b 


_ CO 
>M  A 
DP 
CET 
1Z  M 
CO 

6 I is  K 


75  K=gl,lZR 

iri.23456789 


8 


95 


4 (DN$TY,RA 

ry.1100)  , RA 


-RATIO!  K I * 


TIO,  WALL  ,1 

irioaoQ) , vs 


fiA  L U K)/RH 


I 


J 


0123456789 


I M 1 = 1-1  -- 

11 MK— 1 Z-K 

DE  T { n_=DET(  I MI  )*XMATPXt[I  ZMK  ) 
SUMMON STYt  IZ  ) 

__C(U  85  K=  i * IZ  R 
Kl=  R + l 
_..IZMK=I  Z-K 
SUP=SU M+DNSJY 
PHI  I.  Z ).  = R AT  I 0 ( 

CONT  1NUE 
RETUR'J 
EN  C 


( I ZMK) *CET 
1 Z) ^ SUM 


H,Ph,RHC) 
ALL (100  ) ,P 


01  23456789 


( K L ) 


HI  100  ) ,XMA 


01  23456789 


0123456789 


RXI100) , 


012345678 


_S.UB 

ROUTINE  CU 

r PT  I ( PH_, I H 

). 

CIM 

ENS  ION  PHI 

100  ) 

WRI 

+-* 

rr 

- 

f\j 

s» 

rv 

O 

o 

IP  hi M ) , N = 2 

, IH  ) 

. . 

200  FOP 

MAT  ( T 3 C , P 1 

0.4) 

RET 

URN 

„ 

. 

END 
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1 

1 

1 

1 

2 

2 

2 

2 

2 

3 

3 

3 

3 

2 

4 

4 

4 

4 

4 

"5 

5 

5 

5 

5 

5 

6 

6 

6 

.6 

6 

4> 

6 

7 

7 

7 

7 

7 

7 

7 

e 

8 

8 

G 

G 

8 

6 

8 


pressure  p(y) 

. /m. ^ ) 


Lateral  pressure  p(y) 
(kg. /m. 2 ) 


h 95. 7498 
f36.  C55  9 
>76.  7239 
! 1 1.8083 
43 .362  1 
.71.4397 
>96.0925 
1 7.  3 7 23.  _ 

1 3 5.  32  £4 
>50.0132 
6 1.473  4 
'69.7598 
74.9177' 
76.9954 
>76. C381 
f 72. C 933 
165.2046 
.55.4141 
42  .7695 
27.3125 
’ C9 . 0820 
566.  1250 
>64.4  766 
!28.  1 £36 
► 09 .281 2 
.77  .8086 
' 4 3 • 8 C 8 6 
>07.3125 
;6  6. 3 633 
’26  • 9961 
>83.2461 
l37  • 1 52 3 _ 
.88  .746  1 

9 

565.1367 
3C.  CC39 
? 7 2 . 6 9 53 
13.2422 
.51  .6836 
.88.0430 
122.3555 
>54. 6484 
>64.  9 6 C 9 
>13.3125’ 
139.7  383 
64.2656 
86.9219 
’07.  7383 
!26.  7422 
>43.  957  C 


01  23456789 


1 


9 

9 

9 

9 

9 

9 

9 

10 
loti 
1C 
1C 
10 
J 10 
10 
i ob 
1 c 
1C 
- 1C 
lib 

in 

...  li 

il|2 

113 

nk 

..  ii 
n 

. U|6 
^ 117 
113 

nb 
n 
12 
121 
12 
122 
123 
12  3 
- 12  ^ 
- 12  4 
125 


1 2 3 4 5 6 7 


9059.4102 
9173. 1289 
285.1445 

395.4766 
9504.1523 

611.1953 
71  6.  632  6 
820.4644 
722.7773 

023.5352 
22.7812 

220.5352 
316.  8203 
411. 6641 
50  5.0  761 
1597.0937 
>87  .7266 
777.0000 
864.9297 
951. 543C 

3 6.  6 516 
20. 6626 
1203.6464 
85.1719 
65.4766 
44.5703 

1522.4766 
599.2146 

74. 6CC6 
49 .2539 
22.5859 
94.8164 
965. 9648 
C36.0430 
0 5.  C 6 64 
173.0566 
40.0.273 
05  .9922 
70.9648 
34.9609 
96.  CCCC 
t C . C £ 9 6 

ale[621a2  4 6l8  9 


0123456789 


12 

12 

12 

12 

12 


50  1 .4844 
740.8203 
[7  99.2  656 
656.  832  C 
1913.  5312 
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(3)  Example  No.  3.  Same  as  example  No.  1,  but 


with  alternating  friction,  i.e. 


(° 

y =<  alternating. 

(0.450 


The  results  are  similar  to  example  No.  2.  The • 

effect  of  alternation  between  0 and  0.450  is  similar  to 

y =0.225. 
ave. 

(4)  Example  No.  4.  Simulation  of  the  Reimberts ' 
model.  The  characteristic  functions  are: 

Density- function : 

3 

y = 750  kg/m  (constant) . 

Ratio- function: 


k(y) 


2k 


l+(- 


yk 


R 


)y 


where,- 

k = 0.333 
o 

y = 0.225 
R = 4 . 99  m 


(5)  Example  No.  5.  Simulation  of  dynamic  pressure 
curves.  Same  data  as  in  example  No.  4,  except  for  k(y). 
k (y ) as  before  until  y = 75  m.,  then  decreasing  at  a constant 
rate  of  0.01/m. 


101 


DIGITAL  SIMULATION 
Example  No.  3 


LIM 
...PEA 
1C1  FOR^ 
_CAL 
CA  L 
CAL 
CAL 
CAL 
CAL 
END 


ENS  ION  DNSfT 
C(i  ,101)  1 

1 AT  ( l 3 , P5  • 

1L  SDNS.TYJDN 
L SRATIClK 
L S WA  L L<  WA 
L S PH(  DNST|Y 
3 UTPT1LP 
L EXIT 


Y { 100)  » R A 
H t R HO 

2 ) 

S T V . I H ) 

:1a  T 1 0 , IN) 

L L f I H ) .... 

, RA  T I G , W A 
H , IH  )_ 


ri  0 (1  00)  , WALL  (100  ) ,P 


LL,IH,PH,RHG) 


-i(  100  ) 


SUB 

..DIM 

DO 

25  DNS 
RET 
END 


ROUTINE  SO 
ENS  ION  DNStT 
2 5 J=  1 , 1 H 
r Y1  J ) = 7 5 0 • 
URN 


NSrYl  UNSTY 
Y ( 100  ) 

0 


IH) 


SUB 

ROUTINE”  SR 

A 710 (RAT  10 

,1H) 

...  . 

DI  P 
LU 

ENS  ION  RAT 
25  J = 1 , I H 

I 0 (1  0.0). 

■ - — 

- - - — 

25  RAT 
RET 
END 

1U(  J.)=0  .33 
CRN 

3 

— — — — 

— — ~ — 

SUB 
LI  M 
DU 

25  WAL 
IH  P 
DO 


^UU  TINE  SWA 
ENSIGN  WAL 

25  J=L,IH, 

L (J  ) = C.O 
1=1  H-l 

26  1=2,  IHP1 


26  WAL  L ( I ) = 0 • 4 50 


RET 
EN  C 


UR\1 


LL( WALL, I 
L(  100  ) 

2 


,2 


) 
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SUB 

ROUTINE  SP 

-1(DNSTY,RA 

TI  G,toALL  , I 

H t P H f RHC ) 

CIM 

FNS  ION  DNS 

m 100) , RA 

ri odoo) , w 

ALL (100  ) ,P 

DjETt 

1100|) 

DO 

5 C K=  1 1 I H 

SO  XMA 

IPX  ( K ) = 1 • C 

-R  A T 10  ( K )_* 

WALL! K )/ RH 

0 

j DO 

<5  5 IZi=2,  1H 

1 

4 

5 

01 23456  PlfT 

i 1 ) T 7 J ^ 7 8 9 

01  23456789 

01  23456789 

% i 

01  23456789 

IZM 

1-1  Z-l 

CO 

75  K= I , IZM 

1 

I = K 

f 1 

I M 1 

= 1-  1 

IZM 

| 

1 

!i 

75  CET 

I 1)  -DET(  I M 

1)*XMATRX( 

I ZKK) 

SUM 

= DVI  ST  Y ( IZ  ) 

DO 

85  K=  1 , IZM 

1 

Kl  = 

< + l 

IZM 

|K=i  Z-K 

8 5 SUM 

= SU M+ONST Y 

(1 2MK)*DET 

( K 1 J 

PHI 

IZ  ) =R  A T 1 G( 

11 )*SUM 

: • ; Li 

95  CON 

T IM  LE 

'j 

. i 

■ • -> 

PET 

U KM 

; < 

■ i | • ;l 

END 

SUB 
LI  M 

ROUTINE  CU 
ENS  ION  P HI 

T PT  l(  PH,  IH 
100  ) 

) 

— — *—  f- 

ViRI 

TE( 2,200) 

(PH  M ) , M = 2 

, I H ) 

200  FOR 
RET 
END 

MAT  (T3C,F1 
URN 

0.  5 ) 

- — — — * 

— — — - — — 

( 100  ),XMArftX(100), 

6 7 

1 2345678901  2 3 4 5 6 7 E 
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Lateral  pressure  p(y) 
(kg. /m. 2 ) 


Lateral 

(kg 


1- 


pressure  p(y) 

. /m. 2 ) 


5 

12 
14 
16 
19 
2 1 


4 99  .49976 
734.24976 
83.99951 
C4. 19995 


2 7 
29 
31 


35 

41 

43 

44 

47 

48 
5u 
51 

54 

55 

57 

58 
61 
61 
64 


60.03585 
09.78442 
32.18433 
23  51 .93286 
2531.05444 
8C. 8C298_ 
4 7. C4415 
96.79272. 
3330.54395 
36  00  .29  248 
3 7 1 . 92554 
91.67407 
21. 55C78 
7 1 • 30  C7  6__ 
39.77734 
39  • 52734_„ 
46.9492 2 
96.69922 
93. 39453 
43. 14453 
25.43750 
79.  1 8355. _ 
5 5 .38672 
05.1328  1 
71.54687 
2 1*2  92  97  a 
647  8.21 C54 
67 2 7. 56C54 
6775  .67  187 
7025  .41797 
7054.19531 
7313.94141 
73 44  ,0546  9 
7593.80465 
7615. 50  7 81 
7855 .2578  1 
78/8  .81250 
8126.56250 
.81  34.?].  C94 
8383.56094 
8261.93750 
31.6875C 
22 .22656 
7 1 .97266 
55.29667 


0 1 2 3 4 5 6 


66 

86 

86 

88 


91 
.9  C 

92 
9-2 
95 
95 


230 


1 


0 1 2 3 4 5 6 7 


97 1 
971 
55 

95 

*ib 

*i 

♦ 1 

♦ i 

♦ 1 
_ * i 

♦ 

♦ lb 

♦ 10 

- «i o 

♦ ill 

♦ l 

♦ l 

♦ l 
*1 
♦ltL 

♦ 11 
♦ 11 

♦ ill 

- *1 

♦ 1 
♦ 11 
♦ 12 
♦ 11 
♦ 1|2 
♦ 1 
♦ 1 
+1 
♦ 1|2 
til 


C5.  04 68 7 
81.37109 
31.11715 
6 523  4 
50 .40234 

13  .35  156 
3.  10  156 
5.66406 

165.  41  Cl  6 
15.77734 
169. 527 
0 113.882 
0363.632 
0302 • 160  . 
C551.51Q 
0464. 761 
IP  734.  531 
66  1.921. 

9 1 1.667 
833. 738  _ 
083.488 
1 0 C G . 398 
125C. 146 
1162. 054 

14  1 1.804 
318.85  1 
568.601 
470. 545 
72  C.  695 

1616.466 
1866.216 
76  1.562 
011.312 
900.359 
150. 109 
12  C 34.  56  6 
2284. 736 
2165.574  . 
4 15 .320 

9 


0 1 2 3 4 5 < 


♦1 
♦ 1 
♦ 

♦ 

♦1 

♦1 

♦1 


♦1 
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2541. 588 
2415. 1 Cl 
1 2 6 6 4.  64  7 
24. 269 
2784.019 
? 049.86 3 
2655. 609 
2 761.  580 


DIGITAL  SIMULATION 
Example  No*  4 


DIMENSION  DNST  Y(ICC)  ,RATIO(100)  * WALL ( 100  ) , PH  ( 100  ) 

REAOUtiOl)  IH,RHO __ _ 

ICl  FORMAT!  13  ,F5  .2  f 

CALL  SON  STY  { CNSTY  , IH  ) 

CALL  S RATIO {RAT  10,  IH) 

CALL  S WALL (WALL , IH) 

CALL  S PH (ON^TY, RATIO, WALL,  I hi , P H , R HO  ) 

CALL  OUTPT  11  PH,  IH)  

CALL  £ X I T 

END 


S U BR GU  T I N E ~ SO  NS  T Y { 0 N S T V , I H j 

C IMENS  ION  DNS TY(  IOO)  

DO  25  J = L,  IH 
25  ONSTY(J)= 750.0 
RETURN 
END 


SUBROUTINE'  SRAt  10  { RATIO,  IH) 
DIMENSION  RAT 10(100) 

H Y PE  R=  t 0 . 2 2 5 * C • 33  3 ) / 4 • 99  ' 

DO  2 5 J= 1 , 1 H 

25  RA  T lOTJ  ) -0  .666/  ( U C+H  Y PER* J) 
WRIT  EC  2,100)  {RATICI  J)  ,J=1,IH) 
10 C F 0 RM AT  ( T 20 , F6.4  ) 

RE  TURN 

END*"  ~ " 


SUBROUTINE  S WA  LL(  WA  LL  , IH  ) 

DI-MENS  ION  "WALL  (100) 

DC  25  J= I , I H 
25  WALL(  J ) = 0 . 22  5 
RETURN 

END  
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SUB R CUT  I ME  SPH  ( CN S T Y, R AT  10  » toALL  , I H ,Ph , RH C ) 

T3TMENS  TON  D NSTY  ( ICC),  R ATI  0 (1  00)  , UALL  ( 100  ) ,PH(  100)  ,3(MAT  RXITOOI, 
DETt ) 1100) 

ra so'  k-t  , if 

5C_  XMATRX  (K ) = 1 • O-RAT  10 I K ) *WA L L(  K ) / RH 0 

DO  5512-2 , TH  

DETl  11  = 1.0 


IZ_M1MZ-  1 

CO  75  K = 1 1 1 Z M 1 

I = K+  1 

1 M 1=  I - 1 
IZMK=I  Z-K 

7 5 DE  T ( l ) -DE  T I I Ml ) * XM  A TR  X ( I Z MK ) 
_SUM=DNSTY( 1Z  ) 

~ CO  '85  K = 1 , IZ  P I 

K1  = K +1  

I ZMK=I Z-K 

65  SUN=SUM+DNSTY(  I ZMK)  *DETI  K1  ) 
PH I IZ)  = RATIO(I  Z)*SIM 

95  CONTINUE 


S U6RGU T I NE  CU T P TI  C PH, I F)  " 
DIMENS  ION  PHI  100) 

V'R  I T EC  2,200)  1 P hi  > ) ,>=  2, 1 H ) 
2CC  FORMAT ( T30, F10  .4  ) 

'RETURN 

END 


RETURN 
EN  C 
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(atio  function 
k(y) 

Lateral  pres- 
sure p(v) 
(kg. /m.* ) 

0.655  1 

955.5278 

0.646  6 

1392.3103 

C. 6373 

18C4.3472 

C. 6233 

2193.3657 

0. 6195 

2560. 946C 

0.6  110 

2908.554 2 

0.6027 

3237.5200 

0 .5946 

3549.0813 

0.5  66  7 

3844.3767 

C .5 791 

4124. 4453 

C.571  6 

4390.2656 

0. 5643 

4642.7167 

0 .5572 

4882.6523 

0.5503 

T 5 1 10  .8  20  3 

0.5436 

5327.9375 

0.537  0 

5534.6719 

C. 5306 

5731.6328 

C*  524 3 

551  5.3826 

0.5182 

6056.4646 

0.5122 

6269 .3672 

0.506  3 

6432.5586 

(57500  6 

6588.4609 

0.455  0 

6737.4727 

C.  489  6 _ 

6675. 9561 

C.4842 

7016.3477 

(57479  0 ' _ 

7146.8711 

0 .4739 

727  1.3867 

_ 

7391  .664  1 

0.464  0 

7506.4883 

0.4592 

7616.6055 

0.4545 

7722.2422 

....  .. 

7823. 6445 

0.4  49  5 

7921.0156 

0 .445  3 

80  14.5547 

0.440  9 

8104  731453 

0.436  6 

8190.8672 

0.432  3 

82  73. 9883 

C. 423  1 

8353. 5531 

C.424  C 

84  3 C.  5 16  C 

0*  420  C 

8505. C312 

0.4161 

8576747573 

0 .4  12  2 

864  5 • 1680 

0.4034 

871 1.4609 

0.404  7 

8775.3672 

0.4  Cl  0 

6637. CC35 

C.3575 

8856.4648 

0.3935 

89  532  6 5 54  ’ 

0 . 390  5 

9009.2656 

0 .3870 

9062.7773 

0.383  7 

9114  .4609 

Lateral  pres- 

Ratio  function 

sure  p(y) 

k(y ) 

(kg. /m. 2 ) 

P_L3J°_4 S164.4C  23 


C.3772 

9212. 671 5 

C.3  74  C 

9259.3437 

G.3  70  5 

9304.4766 

0 .3678 

9348 . 1445 

0 .3648 

9390.3867 

0.3618 

9431.2852 

0.358  9 

54  7 C.  8789 

0 . 3 56  0 

9505. 21C5 

C.3  53  1 

9546.3555 

0.3504 

9582.3437 

0.3476 

9617  .2  109 

0 .344  9 

9651.0234 

0.342  2 “ ' 

' 9683.7891 

0.3396 

571 5. 5781 

0.337  C 

5746.41 8C 

0.3345 

9776. 332C 

0.3  32C 

9605.3711 

0 .329  5 

9833  .5586 

ST3I7  1 

9860.9336 

0 .3  24  7 

9887.5117 

0.3224 

5513.3242 

0.3200 

9936. 4C62 

C.3177 

9962.7734 

C.3155 

9986.4727 

0.3132 1 — 

“10009  .496  1 

0.3110 

10031.8945 

0.3  08  9 

10053. 6641 

0.3C67 

1 CC  74.  8554 

G.T04T6'  " — 

10055.4646 

C.3C2  6 

10115. 5156 

crrror  ' — 

10  135.0  352 

0.298  5 

10154  .0273 

0.296  5 

101 72 .5312 

0.294  5 

1 C19C.5547 

07272  6 

1C2C6. C898  " 

C.2S07 

10225.1 756 

T77H8B  

“10241. 6320 

0.2  66  5 

10253.0586 

0 • 28  51“  ~ — 

1027  3 .367  2 

0 .2832 

10285.2891 

0.281 4 

10304.3164 

0.2797 

1C31  6.  572  7 

G.2  77  9 

10  3 3 3.2  656 

0.2  762 

10341.2031 

0.2745 

1 J 360 .7969 

0 .2728 

1 03  7 4 .06b  4 

0.2711 

1C387.0234 

0.2695 

0.267  8 

1 
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DIGITAL  SIMULATION 
Example  No.  5 


01 MENSION  ON  ST Y ( 1 00 ) #R AT  1 01  100 ) , W AL L 1 1 0 C) , PH ! 1 00) 
RE  AD  (1,1  01)  IH,RHO 
101  FOR  RAT  ( I 3,F  5. 2 ) 

CALL  SDNSTY! DNSTY,IH ) 

CALL  SRATIO(RATIO,IH) 

CALL  S WALL  I WALL  , IH) 

CALL  SPH  i ON  STY  f RAT  10  «MA  L L 1 1 H , PH  , RH  C) 

CALL  OUTPTi (PH, IH, RATIO) 

CALL  E XI T 
END 


SUBROUTINE  SDN  STY ( CN STY  , I H ) 
D I ME  NS  1 ON  DNSTY(IOO) 

DO  2 5 J=  1 1 1 H 
25  CNSTYI  J 1 = 750 .0 
RETURN! 

END 


SUBROUTINE  SRA  T I C ( R ATI  0, 1 H) 
CiMENS  ION  RATIO!  100) 

H YPER=  (0.225*0.333  )/4. 99 
DO  25  J=l,75 

25  RATIO!  J)  = C.  666 / ( 1 . O+HYP  ER  * J) 

DO  35  J=  76 , 1 H 

35  RATIO!  J)=RATIC<  J-l)  -0.  Cl 

WRITE!  2,100)  (RATIO!  J)  ,J=1,IH) 
100  FORMAT ( T20 , F6 .4 ) 

RETURN 

END 


SUBROUTINE  S WA  L L(  WA  L L , IH  ) 
DIMENSION  WALL ( 100  ) 

DO  25  J=1,IH 
25  WA  LL(J ) = C . 2 2 5 
RETURN 
END 
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SUBROUTINE  SPH  I CN  S TY,R AT  10  » WALL  ,IH ,PH,RHC) 

DIMENS  ION  DNSTY(100)fRATIO(100>  ,WALL(100),PH(100),XMATRX(100), 

DETt  1100) 

CO  50  K=  1 f I H 

50  >M  ATRX  IK ) = 1 «0-RAT  10  ( K)  *WALL(  K)  / RHO 


DO  95  12=2, IH 
DET  I 1)  = 1.0 
IZM  1 = 1 2-  1 
CO  75  K=  1 » IZ  M 1 
l=K+l 
I M 1=1-1 
IZ  PK=I  Z-K 

75  DETI  I)  =DET  1IMI)*XMATRXI  IZMK ) 
SUM=DNST Y I I Z ) 

CO  85  K=  1 , IZ  P 1 
Kl=K+i 
I ZMK=I  Z-K 

85  SUK=SUPH-DNSTY(I  ZMK  ) *DET{  K1  ) 
PHI  I Z) =RAT 10(1 Z)  * SUM 

95  CONTINUE 
RETURN 
END 


SUBROUTINE  CUT PT1 { PH , I H, R AT  10 ) 

DIMENS  ION  PHI  100) , RATI  01100)  ,21100) 

ViR  I T E(  2 , 200  ) I PHI  P),P=2,  IH) 

200  FORMAT  I T3  0, F10  *4 ) 

Zi  1)  = 1.0 
DO  5C  1=1, IH 
50  ZI  1*1)  =Z 1 1 ) + l.  C 

CALL  MAXI  PHI  2) ,98,VALMA>,  I SUB) 

CALL  1 INIPHI  2),98,VALPIN,ISUB) 

CALL  PL0T(0,98,2,VALMAX  ,VA  LPI  N,  Z,  PH,  R * TI  0) 

RETURN 

END 
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Ratio  function 
k(y ) 

0*656  1 
0 • 646  6 
0*6373 
0*6283 
0.6195 
0*6110 
0.6027 
0.5946 
0.5  86  7 
0.5791 
0.571  6 
0.5643 
0.5572 
0.550  3 
0.543  6 
0 . 537  0 
0.5306 
0.5243 
0.5182 
0.512  2 
0.506  3 
0.5006 
0.495  0 

0.4696 
0.4842 
0.4790 
0.4739 
0.468  9 
0.464  0 
0.4592 
0.4545 
0.4499 
0.4453 
0.4409 
0 .436  6 
0.432  3 
0.4281 
C.424C 
0. 420  C 
0.416  1 
0 .4  12  2 
0.408  4 
0.4C47 
0.401  0 
0.3575 
0.353  5 
0 .350  5 
0.387  0 
0 .383  7 


Lateral  pres- 
sure p(y) 
(kg./m.^ ) 


955.5278 
1392.3103 
1804.3472 
2153.3657 
2560.  54  8 C 
2908. 5542 
3237.5200 
3549.0813 
3844.3767 

4124.4453 
435C.2656 
4642. 7187 
4882.6523 
5110.8203 
5327.9375 
5534.6719 
5731.6328 
5515.3  828 
6058.4646 

6269.3672 
6432.5586 
6588  .4609 
6737.4727 
6879.9961 
7016.3477 
7146.8711 
7271. 8667 
7391.6641 
7506  .4883 
7616.6055 
7722.2422 
7823.6445 
7921.0156 
80  14.  5547 

8104.4453 
8190.8672 
8273.9883 
8353.9531 
843  C. 9180 
8505.0312 
8576.4023 
8645.1680 

871 1.4609 

8775.3672 
683  7.  CC35 
8896.4648 
8953.8554 
9009.2656 
9062.7773 


Ratio  functi 
k(y) 


0.3804 
0.3772 
0.374C 
0.3705 
0.3678 
0.3648 
0.3618 
0.358  9 
0 . 356  0 
0.3  53  1 
0.3504 
0.3476 
0.3449 
0.342  2 
0.3396 
0.3370 
C.  334  5 
0.332  C 
0.3295 
0 .327  1 
0 . 324  7 
0.322  4 
0.320  0 
0.3177 
0.3155 
0 . 3 13  2 
0 .303  2 
0.293  2 
0.2832 
0.2  73  2 
0.2632 
C.2  53  2 
0.2432 

0.233  2 
0.223  2 
0.2132 
C.2C32 
C.  193  2 
0.1632 
0.1732 
0.163  2 
0.153  2 
0.143  2 
0.133  2 
C. 123  2 
0.113  2 
0.1032 
0 .0932 
0.083  2 
0.0732 


Lateral  pres 
sure  p(y) 
(kg. /m.* ) 


9114.4609 
9164.4023 
9212.671  5 

9259.3437 
9304.4766 
9348.1445 
9390.3867 
9431.2852 
947C.8789 
9 509.  2 1 C 5 
9546.3555 

9582.3437 
96  17.2109 
9651.0234 
9683.7891 
9715.5781 
9746.41 6C 
9776. 332C 

9805.3711 
9833.5586 
9860.9336 
9887.5117 
9913.3242 
9936. 4C62 
9962.7734 
9735.9375 
9506.0703 
9272.9297 
9036.2266 
8795. 707  C 
8551. C7C3 
8302.0547 

8048.3711 
7789.7266 
7525.8242 
7256.3672 
6961. C46  9 
6699. 5391 
6411.5312 
6116.6836 
5814.6602 
5505.0977 
5187.6367 
4661. 9141 
4527.5234 
4184.0781 
3831.1609 
3468.339  1 
3095.1707 
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(6)  Example  No.  6.  Simulation  of  dynamic  pressure 


es.  Same  as  example  No.  5,  but  k (y)  = 0.100  (constant) 


for  y > 75  m. 
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DIGITAL  SIMULATION 
Example  No*  6 


DIMENS  ION  DNSTY { IOC) , RATI  0 (100) , WALL (100)  ,PH(  100 ) 
READd  ,101)  IH  ,RHO 
1C1  FORMAT (I  3, F5.2 ) 

CALL  SDN  STY ( DN  STY , I H ) 

CALL  S RAT  I 0 ( RA  TI 0, IH) 

CALL  S WALL(WALL,IH) 

CALL  SPH( DNSTY , RATI C , W ALL , I H, PH, RHO ) 

CALL  3UTPT1(PH,IH, RATIO 

CALL  EXIT 

END 


SLBROUT INE  SDNSTY ( DNST  Y, I H) 
DIMENSION  DN  STY  ( 100  ) 

DO  25  J=  1 , I H 
25  DNSTYI  J ) = 750.  0 
RETURN 
END 


SUBROUTINE  SRAT10  (R  AT  10  , I H) 
DIMENSION  RATI  0 ( l 00) 

FYPER= ( 0 • 225*0. 333) /4. 99 
CO  25  J= 1 ,75 

25  RATIO!  J)  = 0*666/(  1.0+HYPER*J) 

DO  35  J=  76 , IH 
35  RATI  Cl  J)=C.100 

WRITE!  2, 1 CC)  (RATIO (J) ,J=1 ,IH) 
100  FORMAT  (T20,F6*4) 

RETURN 

END 


SUBROUTINE  SW ALL ( W ALL , I H) 
DIMENS  ION  WALL(  100) 

CO  25  J ~ 1 , 1 H 
25  WALL  (J  )=  0.225 
RE  TURN 
END 
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SUBROUTINE  SP H( DNSTY  ,RATI C, WA LL ,1  H , PH, RHO) 

DIMENSION  DNSTY (100)  ,RAT  1 01  100 ) , WALL  1 1 00)  , PH ( I 00) , XMA TRX l 100 ) , 
DET  ( 1100) 

DO  5C  K=  1 » l H 

50  XMA TRX ( K ) = 1. 0— RATI  01K)*WALL(K)/ RHO 


CO  95  12=2, IH 
DETI1)  =1.0 
12  Ml=l  2-1 
DO  75  K=1 , I 2M1 
I = K + 1 
1M  1=1-  1 
I2MK=12-K 

75  DE  T ( I)  = DET  ( IM  1 ) *XMA  TR  X! I2MK) 
SUM=UN STY  1 1 2) 

00  85  K= l , I 2MI 
K l=K+l 
I2MK=I2- K 

85  SUM= SUM +DNSTY!  I2MK  )*0E TIKI ) 
PH  II  2)  =R AT  10  ( 12  )*SUM 

55  CONTINUE 
RETURN 
END 


SUBROUTINE  0 LTPT  1 ( PH  , I H ,RAT  I C ) 
DIMENSION  PH (100) .RATIO! 100 ) ,21 100 ) 
WRITE!  2 » 2 C C ) { PH  (M)  , H=2  , I H) 

200  FORMAT!  T30.E10.  A) 

2 ( 1 ) =1  .0 
CO  50  1=1, IH 
5C  21 1 + l)  =2  ( I ) + 1.0 

CALL  MAX! PH  12) ,98,VALMAX,ISUB) 

CALL  MIN! PH(2),98,  VALMIN.ISUB) 

CALL  P LOT (0,98, 1,VALMAX, VALMIN,2,PF) 

RETURN 

ENC 
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Ratio  function 
k(y) 


Lateral  pres- 
sure p(y) 
(kg./m.S) 


0.6561 

955. 5278 

0 .6466 

1392.3103 

0.6  37  3 

1804.3472 

0.6283 

2193.3657 

0.6195 

2560.9480 

0.611  0 

29C8.5542 

0.6C2  7 

3237. 52CC 

0.  5946 

3549.0813 

0.586  7 

3844.3767 

0 .5  79  1 

4124.4453 

0.571  6 

4390.2656 

0.5643 

4642.7187 

0.5572 

4882. 6523 

0.5503 

5110. 82C3 

0.543  6 

5327.9375 

0 . 5 37  0 

5534.6719 

0.5  30  6 

5731.6328 

0.524  3 

5919.3828 

0.5182 

6098.4648 

0 . 5 12  2 

6269.3672 

0.5C63 

6432.5586 

0.5C06 

6588.4609 

0 .495  0 

6737.4727 

0.4896 

6879.996  1 

0.4842 

7016.3477 

0.479C 

7146.8711 

0.4739 

7271.  8867 

0.4689 

7391. 6641 

0.4640 

7506.4863 

0.459  2 

7616.6055 

0.4545 

7722.2422 

0.4499 

7823.6445 

0.445  3 

7921.0156 

0.4409 

8014.5547 

0 . 4 36  6 

8 10  A. 4453 

0.4  32  3 

8 190. 8672 

0 .4  28  1 

8273.9883 

0.4240 

8353.953  1 

C. 420  0 

8430.9180 

0.416  1 

85C5.0312 

0.4122 

8576.4023 

0.4084 

8645. 168C 

0.4047 

8711.4609 

0.4010 

8775 .3672 

0.3975 

8837.0039 

0.393  9 

8896.4648 

C. 390  5 

8953. 8594 

0.387C 

9009.2656 

0 . 3 £3  7 

9062.7773 

0.380  4 

91  14.4609 

Ratio  function 
k(y) 

Lateral  pres- 
sure p(y) 
(kg. /m.* ) 

0.3772 

9164.4023 

0 . 3 74  0 

9212.6719 

0.3709 

9259.3437 

0.367  8 

9304.4766 

0.364  8 

9348.1445 

0.361  8 

9390. 3£67 

0.3589 

9431.2852 

0.356  0 

9470.8789 

0.353  i 

9509.2109 

0.3504 

9546.3555 

0.3476 

9582.3437 

0.3449 

9617.2109 

0.342  2 

9651.0234 

0.3396 

9683.7891 

0.3370 

9715.578  1 

0.334  5 

9746.4180 

0.332  0 

9776.3320 

0.3295 

9805.3711 

0.3271 

9833.5586 

0.324  7 

9860.9336 

0 . 322  4 

9887.5117 

0.32)0 

9913.3242 

0.3177 

9938.4062 

0.3155 

9962. 7734 

0.3132 

3210. 5771 

C.1C0C 

3271.0999 

0. 100C 

333  1.3499 

0 . 100  0 

3391.3284 

0.1000 

3451.0356 

0.1000 

3510.4775 

C.1C0  0 

3569.6450 

C.  1000 

3628. 5513 

0. 100C 

3687. 1899 

0 . 100  0 

3745.5623 

0.100  0 

3803.6719 

0.1000 

3861.5215 

0.1  coo 

3919.1113* 

C.1C0  0 

3976.4387 

C.  1000 

4033. 5C93 

0. 100C 

4090.3198. 

0. 100  0 

4146.8789 

0.100  0 

4203.1758 

0 . 1 00  0 

4259.2266 

0.1 00  0 

431 5. 0195 

0.  1C0  0 

43  7 C.  5625 

C.  ICOC 

4425. 8555 

0. 1C0C 

4480. 9C23 

0. 100  0 

4535.6953 

0.100  0 

11)  i. 


CHAPTER  6 


MODELS  OF  HIGHER  COMPLEXITIES 

6.1  FORMAL  GENERALIZATIONS  OF  THE  ONE -DIMENSIONAL 

MODELS 

The  one-dimensional  models  can  be  generalized 
formally  into  higher-dimensional  models.  The  derivations 
and  the  results  become  more  complicated  and  therefore  lose 
some  of  their  practical  value.  A more  practical  way  of 
generalization  is  presented  in  section  6.2.  In  the 
following  formal  generalizations  to  two-dimensions  are 
presented. 


6.1.1  Analytic  Model 

In  the  two-dimensional  problem  all  the  variables 
are  considered  functions  of  either  the  depth  coordinate 
y alone,  or  both  y and  the  transverse  coordinate  x. 

By  a similar  development  to  that  presented  in 
section  3.1.2,  a solution  parallel  to  Eq.  (3.13)  is  obtained 
in  the  form 


p <p 


'y>  - \f 

o 


{k (x,y ) [exp ( 


-j S (x,y ) dy ) [Jy  (x,y)  exp  {Jfc 
-Jy  (x,y)  exp  (x,y)  dy)  dy  |y=Q ]\}dx 


y)<ay)  [/y  (x , y ) exp  (/ 6 (x , y ) dy ) dy  - 

(6.1) 
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where 


(6.2) 


6.1.2.  Algebraic  Model 

Similarly  to  the  development  in  section  3.2.1, 
a net  is  formed  over  one-half  of  the  rectangular  vertical 
section  of  the  infinite  bin.  Starting  from  the  axis  of 
symmetry,  v subdivisions,  each  of  width  Ax,  are  made  along 
the  x-direction;  the  vertical  partition  remains  the  same 
as  in  the  one-dimensional  case.  Each  smaller  rectangle 
belonging  to  the  net  is  double-indexed  by  i and  j where 
l£i<u,  1<  j <v , and  so  are  the  corresponding  values  of  the 
local  parameters.  Thus,  the  parameters  are  not  represented 
as  continuous  functions,  but  as  arrays  of  discrete  points 
which,  in  turn,  can  be  represented  in  the  form  of  v-by-u 
matrices.  Then  the  two-dimensional  solution  comes  out  in 
the  form 


v 


n 


p = p K . y , .det  Sub . [I- ( AK)  . ] 
nv  jLj  nj  Z-/  n+l-i,j  l j 


j=i  i=i 

or,  in  terms  of  k . 

n: 


v n 


n 


where 


0 


0 


A ~ .k  ~ . 
n-2 ,j  n-2 , j 


(Ak) 


0 


6.1.3  The  Characteristic  Functions 

(1)  The  density-function.  The  behavior  of  the 
density-function  along  the  y-direction  was  discussed  in 
section  4.2.  The  distribution  of  the  density  along  the 
x-direction  may  be  assumed  to  be  uniform  for  most  practical 
cases.  However — to  be  exact--it  is  probable  that  the 
density  is  larger  towards  the  center  line  of  the  bin.  If 
the  density  distribution  in  the  x-direction  is  assumed  to 
behave  parabolically , and  the  density  distribution  in  the 
y-direction  is  according  to  Eq.  (4.2),  then  the  resulting 
two-dimensional  density-f unction  is  illustrated  in  Figure  6.1. 


function  is  too  complicated  and  is  not  necessary.  The 
geometric  illustration  is  qualitative.  A more  precise  diagram 
can  be  constructed  in  any  particular  case  and  the  entries  of 
the  density-matrix  can  be  measured  directly  from  the  diagram 


The  derivation  of  an  analytic  expression  for  this 
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and  used  when  evaluating  Eq.  (4.2). 


(2)  The  ratio-function.  The  behavior  of  the 

ratio-function  along  the  y-direction  was  discussed  in 
section  4.3.  The  behavior  along  the  x-direction  is  even 
more  complicated.  There  is  no  experimental  data,  nor 
theoretical  knowledge  of  this  factor.  If  assumed  constant 
with  respect  to, the  x-direction,  a typical  two-dimensional 
ratio-function  is  illustrated  in  Figure  6.2;  otherwise, 
it  may  become  as  complicated  as  Figure  6.3  illustrates. 

(3)  The  friction-function.  The  friction-f unction 
remains  basically  one-dimensional,  a case  treated  in 
section  4.4. 
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->T(x,y) 


Figure  6.1. — Two-dimensional 
function  over  a vertical  s 
of  the  infinite  bin. 


k(x,y) 


Figure  6.2. — Two-dimensional  hyperbolic 
ratio  function  over  a vertical  section 
of  the  infinite  bin. 


k(x,y) 


Figure  6.3. — Two-dimensional  dynamic 
ratio  function  over  a vertical  section 
of  the  infinite  bin. 


6.2 


COMPUTER-ORIENTED  GENERALIZATIONS 


As  seen  in  the  previous  section,  formal  generalizations 
of  the  one-dimensional  models  lead  to  somewhat  complex 
expressions  difficult  to  handle  in  practice.  An  alternate 
approach  is  to  keep  the  basic  models  one- dimensional  and 
reduce  higher-dimensional  problems  to  series  of  one- 
dimensional problems.  It  is  best  to  illustrate  this  technique 
by  a concrete  example. 

6.2.1  Example:  A Square  Cross-Section 

(1)  Janssen's  conditions.  Consider  a bin  of  a 
square  cross  section  (Fig.  6.4k.a)with  sides  of  length  a. 

Assume  static  (Janssen's)  conditions.  All  parameters  are 
uniform  and  constant.  For  reasons  of  symmetry  treat  one 
quarter  only,  and  within  the  quarter  treat  one-half  only. 

Fig.  6.4.b.  Subdivide  it  into  a few  sections  of  width  Ago 
each.  The  length  of  each,  p^,  vary.  Each  section  can  be 
treated  as  a single  section  of  an  infinite  bin  with  p=p^, 
and  the  one-dimensional  model  can  be  used  repeatedly  to 
calculate  the  pressure  section  by  section.  For  example, 
if  a = 16',  Ago  = 1',  this  method  requires  eight  repeated 
calculations  with  pi  ,P2 , . . . , pe / compared  to  just  one 
required  by  the  conventional  practice.  However,  the 
resolution  obtained  is  quite  interesting.  Compare  the 
results  of  the  two  techniques. 
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A conventional  calculation  would  give  a rectangular 
pressure  diagram,  Figure  6. 5. a,  i.e.  uniform  lateral  pressure 
across  the  wall.  Examine  the  results  of  the  new  method. 

The  lateral  pressure  for  each  section  is  given  by 

P • Y yk 

Pi  = -j]-  [l-exp{-  (— ) y } ] . 

Only  varies  linearly  from  section  to  section.  It  can  be 

seen  that  p.  is  basically  directly  proportioned  to  p. 

piy 

because  of  the  dominating  — jj—  term.  This  would  give  a 
triangular  pressure  profile.  However,  this  is  modified 
slightly  by  inverse  proportionality  due  to  the  appearance 
of  p^  in  the  denominator  of  the  negative  power  of  e.  The 
combined  effect  is  a truncated  apex  of  the  triangle,  and 
the  profile  is  closer  to  a parabola.  Fig.  6.5*b. 

Now  it  can  be  seen  that  the  old  method  simply 
averages  the  triangular  profile  of  the  new  method.  This 
"averaging"  ignores  the  fact  that  at  the  middle  the  pressure 
may  be  nearly  twice  as  high  as  calculated  by  the  old  method. 

A practical  suggestion  is  to  use  some  compromise,  perhaps  as 
in  Fig . 6 . 5. c. 

Surcharge  can  be  added  in  a similar  way  by  segmenting 
the  pyramid  over  the  square  cross-section  into  sections 
corresponding  to  the  sections  beneath.  The  method  can  be 
extended  easily  to  cross-sections  of  arbitrary  shapes. 
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is  uniquely  associ-ated  with  two  integers  i,  j,  where 
1 £ i _<  p,  1 £ j £ q and  pq  is  the  cardinal  number 
of  the  set  A denoted  by  card  A.  An  example  of  such  a 
net  would  be  a net  of  uniform  partitions  in  both 
directions  consisting  of  intervals  of  length  d,  where  d 
is  the  diameter  of  the  spherical  idealized  kernels. 

This  net  is  adopted  from  now  on. 

(2)  The  image  set.  The  associated  characteristic 
half-pile  is  a right  triangle  T.  It  is  to  be  noted  that 
certain  regions  of  the  triangle  and  rectangle  overlap. 
Call  the  collection  of  the  images  of  all  a e A (included 
within  the  boundaries  of  T)  the  image  set  B.  The 
boundaries  of  T can  be  established  from  two  assumptions: 
(i)  The  equality  of  areas  bounded  by  R and 
T 4/  (i.e.,  card  A = card  B) . 

(ii)  The  imaginary  angle  of  repose  equals  the 
natural  angle  of  repose. 

The  height  h and  the  base  r of  T are  therefore 


h = V2Hptan  0 

(7.1) 

r = ~\/2Hpcot  0 

(7.2) 

It/ An  area-preserving  mapping,  implied  by  the  assumption 
of  a constant  density. 


132 


(2)  Other  conditions.  If  the  various  parameters 


are  not  uniform,  the  calculation  of  each  section  simply 
utilizes  the  local  parameters  independently  of  neighboring 
parameters.  The  total  pressure  diagram  is  put  together 
after  all  sections  are  calculated.  The  technique  of 
Appendix  B can  be  used  in  this  connection  to  give  three- 
dimensional  representations  of  lateral  pressure  diagrams 
over  the  walls. 

6.3  FEEDBACK  AND  COUPLING  MECHANISMS 

6.3.1  Introduction 

Many  of  the  factors  entered  in  the  one-dimensional 
models  are  interrelated  through  various  direct  and  indirect 
relationships.  For  example,  the  density,  y,  is  a function  of 
many  factors,  including  the  vertical  pressure  q,  i.e. 

y = y (q)  • 

However,  q is  or  course  always  a function  of  y,  i.e. 

q = q (y)  • 

Both  are  functions  of  the  depth  variable  y,  i.e. 

y = Y(y,q) 

and 

q = q (y  fY  (y  *q) ) • 

I * 
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Examining  the  last  function,  it  can  be  seen  that  q is  ulti- 
mately dependent  on  itself  through  an  indefinite  loop  that  may 


be  written  symbolically  as, 

q = q (y  rY  (y  (y /Y  (y  rq  ( )• 

In  diagrammatic  form  the  relationship  between  q,  y,  and  y can  be 
presented  as  follows, 


where  an  arrow  pointing  from  a box  to  another  indicates  the 
dependence  of  the  variable  in  the  first  box  on  the  variable 
in  the  second.  Here,  q and  y are  dependent  variables  and  y 
is  an  independent  variable. 

The  following  conventions  are  adopted: 

Dependent  variable  - at  least  one  arrow  originates 

from  box. 

Independent  variable  - no  arrow  originates  from  box. 
(Box  only  receives  arrows). 

6.3.2  Modified  Incidence  Matrix 

When  many  variables  are  involved,  larger  and  more 
complex  diagrams  can  be  constructed.  In  a given  grain  pressure 
problem,  let  N be  the  number  of  dependent  variables, 
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VN,  and  M the  number  of  independent  variables,  , V2 , . ..,  V^. 
Various  relationships  among  variables  (both  dependent  and  indepen- 
dent) may  hold.  For  example, 


V.  + V. 
i D 


V.  f V. 
i 1 


v.  «-  V. 
1 1 


V.  + V, 
l k 


V.  / V, 
l ' k 


for  some  or  all  j ^ i 


for  some  or  all  j 


for  some  or  all  j / i 


for  some  or  all  k 


for  some  or  all  k. 


In  addition,  transitivity  holds, i.e.,  and 

imply  -*■  . 

In  order  to  organize  a problem  it  is  best  to  represent 
the  interrelationships  of  variables  in  a modified  incidence 
matrix  representing  the  directed  graph  associated  with  the 
appropriate  box  diagram.  The  ' s are  listed  horizontally  and 
vertically,  followed  by  the  ' s . The  numbers  0 through  3 are 
entered  in  the  intersecting  squares  according  to  the  given 
relationships  between  the  variables.  This  is  best  illustrated 
by  an  example. 
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Example 


Given:  N = 5 , M = 2. 


Relationships:  V^r  V^*«-  V4  r + V5 , -*  V2 


V2  + V3 ' V3"  V2 ' V5  + V1 


The  modified  incidence  matrix  is : 


V1 

V2 

V3 

V4 

V5 

V1 

0 

0 

0 

2 

2 

V2 

1 

0 

0 

0 

0 

V3 

0 

1 

0 

0 

0 

V4 

2 

0 

0 

0 

0 

V5 

2 

0 

0 

0 

0 

V1 

0 

0 

0 

0 

3 

V2 

3 

3 

3 

3 

3 

The  usefulness  of  this  matrix  is  in  that  it  codes 
the  interrelationships  of  the  entire  system  in  a simple  numeric 
form  that  can  be  easily  fed  to  a digital  computer,  analyzed 
further  and  help  in  the  organization  of  the  calculations  of 
certain  coupling  effects  described  in  the  next  section. 
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6.4  AN  EXAMPLE  OF  A COMPLEX  COUPLING  MECHANISM 

The  basic  parameters  of  a grain  pressure  problem 
are  the  density,  y;  the  vertical  and  the  lateral  pressures, 
q and  pi  and  the  frictional  stress,  s — q,  p,  and  s are 
interrelated  via  the  ratio  and  the  friction  functions,  k and  X. 
The  interdependence  of  these  parameters  can  be  represented 
in  a box  diagram  as  follows. 


Using  this  box  diagram,  all  the  possible  coupling 
effects  due  to  changes  (increase  or  decrease)  in  the  value 
of  each  parameter  are  enumerated  as  follows: 
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+ 

Y = 


==>  q 


+ 

q = 


==>  Y 


+ 

s = 


+ 

P = 


===>  s 


+ 

q = 


q ===>  y 


s ===>  q 


p ===>  s 


q ===>  p 


Now,  consider  the  "chain  reaction"  (fig.  6.6)  resulting  from 
a single  instantaneous  change  in  A,  a likely  occurrence  under 
dynamic  conditions.  Suppose  that  X increases,  which  causes 


Figure  6.6  A path  of  a "chain  reaction"  in  a grain  pressure 
system 


enter  here 


/ + ^ + 

• y ===>  q 


\ 


v 

\ 

+ ^ . 

q > y 


s+  ===■ 


+ ====>s  s+ 


\ V 

\ N 


q+  ===>X p+ 


y ===>  q 


\ 


\ 


\ 

\ _ 

q ===>  Y 


\ 


s_  ===>~~q+ 

\ 

\ 

p”  ===>  s 

V 

\ 

q-  ===>X p“  / 


\ 


/ 
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s to  increase.  Now  enter  figure  6.6  at  s+  and  follow  all  the 
resulting  changes  by  drawing  connecting  arrows.  As  can  be 
observed,  the  entire  system  is  interconnected  by  this 
chain  reaction  and  finally  loops  to  the  starting  point,  to 
start  a new  round.  Computationally  this  means  that  the 
entire  system  has  to  be  recomputed  iteratively  until 
relaxation  is  obtained. 

Suppose  now  that  X decreases,  which  causes  s to  de- 
crease. Here  the  figure  should  be  entered  at  s“,  If  and 
X alternate,  as  may  occur  under  dynamic  conditions,  the  system 
oscillates  between  two  entry  points  s+  and  s , while  two  disparate 
chain  reactions  continue  to  propagate  throughout  the  system. 

This  analysis  shows  that  a grain  pressure  system  is 
indeed  inherently  unstable.  Even  under  static  conditions,  small 
changes  in  X for  example,  due  to  external  effects,  may  cause 
unpredictable  changes  throughout  the  entire  system.  This  may 
account  for  some  of  the  observations  of  small  changes  in 
pressures  under  static  conditions. 
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CHAPTER  7 


TOPOLOGICAL  MODELS 


7 . 1 THEORY 

7.1.1  Introduction:  The  Associated  Characteristic  Pile 

The  theory  of  grain-pile  transformations-  (GPT)  is 

based  on  pure  mathematical  ideas,  as  well  as  on  the 
empirical  properties  of  grain  under  the  influence  of 
gravity.  It  is  well  known  that  a mass  of  grain  supports 
itself  in  certain  unique  shapes (such  as  a cone,  a pyramid, 
etc.)  having  constant  slopes,  under  natural  conditions. 

The  configuration  of  grain  in  a deep  bin  may  be  considered 
as  a geometrical  distortion  of  the  "natural  formation" 
(characteristic  pile)  which  the  same  mass  would  have  taken, 
if  not  constrained  by  the  walls.  Thus,  a stored  mass  of 
grain  in  a deep  bin  has  a potential  tendency  to  transform 
into  a predictable  geometric  figure.  The  pressure  exerted 
on  the  bin  walls  may  be  regarded  as  a result  of  this  intrinsic 
tendency. 

Therefore,  the  geometry  of  a given  bin  together  with 
the  properties  of  the  stored  grain  are  uniquely  associated 
with  some  characteristic  pile.  For  example,  a cylindrical 
bin  is  associated  with  a right  circular  cone,  the  dimen- 
sions of  which  are  dependent  on  the  dimensions  of  the 
circular  cylinder,  and  its  slope  is  dependent  on  the  intrinsic 
properties  of  the  grain.  Factors  such  as  the  angle  of  internal 
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where 


H - height  of  R 

p - width  of  R 

0 - natural  angle  of  repose. 

Thus  T is  determined  in  terms  of  the  bin  geometry 
and  the  grain  properties.  Extend  the  net  to  include  T too, 
such  that 

(i)  the  lower  right-hand  side  vertex  is 
double-indexed  by  1,1 

(ii)  the  lower  left-hand  side  vertex  is  double- 
indexed  by  r',  1 

(iii)  the  upper  vertex  is  double -indexed  by  l,h* 

where 

r'  = g-  and  h*  = jj  . 

(3)  Definitions . All  terms  used  in  this  section 
are  to  be  understood  precisely  as  defined  below,  and  not  as 
interpreted  in  general  usage. 

a.  Let  a^jSA.  A neighborhood  of  a^j,  designated 
by  N[a^j],  is  a non-empty  subset  X of  A,  such 
that 

(i)  a . . eX 
13 

(ii)  there  exists  an  x eX3 x ^a. . 

mn  irur  lj 

(iii)  x eX=>x'  eX,  where 
mn  T uv 

| i - u | <J  i - m | 

| j - v|<|  j - n | . 
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b. 


A neighbor  point  of  a point  a^eA  is  a 

point  x j^a.  .eN[a.  .]  such  that 
mn  13  lj 

| i - m | =1  or  0 

| j - n | =1  or  0 

c.  An  interior  point  of  A is  a point  a^eA  that 
has  at  least  four  neighbors. 

d.  A boundary  point  of  A is  a point  a^eA  that 
has  less  than  four  neighbors. 


7.1.3  Grain-Pile  Transformations  (GPT) 

Since  card  A = card  B,  there  exists  a non-empty 
finite  set  ITl  (card  111  = (pq) I ) such  that  if  Me  m then 

(i)  M is  a mapping,  M:A+B 

(ii)  M is  one-to-one 

(iii)  M is  onto. 

(Note:  A mapping  M from  a set  A to  a set  B is 

to  be  one-to-one  and  onto,  if  for  every  element  belonging  to 
A there  exists  a unique  element  belonging  to  B,  and  vice  versa.) 

Among  all  Me  m,  seek  a unique  mapping  Mq  such  that 
particular  physical  constraints  are  satisfied.  The  physical 
constraints  are  presented  in  the  form  of  propositions  as  follows: 
PROPOSITION  1:  The  physical  ordering  of  the  initial 

set  A over  R is  preserved  in  some  unique  manner  (to  be  dis- 
cussed later,  see  section  7.1.4)  in  the  image  set  B over  T. 
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PROPOSITION  2:  Given  any  pair  (a  — / ^o^aij^eMo' 

it  is  possible  to  find  a continuous  curve  (of  non- 

negative length)  joining  a^^  and  its  image  point  M^a^j), 
card  (^  = pq,  such  that  if  all  the  a^j's  physically 
travel  along  the  corresponding  C^'s  when  a physical 
transformation  actually  occurs  6/,  some  physical  factor  6/ 
is  a minimum  with  respect  to  any  other  set  0'  (card  Q'  = pq) 

of  continuous  curves  C! . each  joining  a. . and  M (a. .). 

13  J lj  o ij 


7.1.4  Some  Intermediate  Results 

(1)  The  decreasing  character  of  the  dynamic  ratio- 
function.  In  line  with  proposition  1,  assume  the  following 
properties  of  the  mapping  Mq: 


5/  Such  a physical  transformation  never  does  occur  in  reality 
during  the  life-time  of  the  structure.  However,  when  the  walls 
yield  to  some  sufficient  extent,  the  transformation  is  "nearly" 
started.  The  forces  acting  on  each  kernel  at  this  instant 
are  tangent  to  the  C. . curves  at  the  initial  points  a. .. 

From  the  point  of  vie^  of  these  forces,  it  is  an  1-! 

immaterial  coincidence  that  at  the  successive  instant  the 
walls--rather  than  proceeding  to  yield--tend  to  return.  This 
situation  is  typical  to  the  dynamic  effects  during  discharge. 

6_/  In  example,  work  or  time.  The  latter  choice  suggests  a 
notable  resemblance  to  the  classical  brachistochrone  problem. 
This  one,  however,  seems  to  be  much  more  complicated  and  may 
be  termed  therefore  as  a "finite  multi-brachistochrone" 
problem. 
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( i ) Given  a . . e A , then  w N [M  ( a . . ) ] : 
ij  v o 1 j 

(ia)  if  a^j  is  a boundary  point,  a a point  t 


-1 


in  every  N[MQ(a^j)]  3 Mq  A(t)  is  a neighbor  of 


and 


ai  j * 

(ib)  if  a^j  is  an  interior  point,  a points 
t^  and  t2  in  every  N[MQ(a^j)]  3 Mq  ^(t^) 

M ^(t0)  are  distinct  neighbors  of  a. .. 
a a non-empty  proper  subset  ICM^  (a^j, 

Vaij))eI  =**  aij  = Mo(aij>- 

(iii)  Boundary  conditions  for  the  image  set: 


(ii) 


(iiia) 

(au- 

Van>)eI 

(iiib) 

(alh'  ' 

Vaih’))eI 

(iiic) 

M (a 

) =»  b , , . 

o'  pq 

r 1 1 

If  (i) 

is  not 

satisfied  for  some 
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then  such  an  image  point  is  called  a point  of 
discontinuity,  or  a singular  point  of  B. 

These  properties  almost  determine  the  mapping  Mo« 

The  second  proposition  will  be  used  to  complete  the  determination 
of  Mq.  In  light  of  proposition  2,  ^ is  a family  of  mutually 
non-intersecting  catenary- type  curves. 

Knowing  and  Q , it  is  possible  to  find  the 
tangent  of  the  angle  (a..^)  between  the  tangent  line  to 
at  a. . and  the  vertical  direction.  The  ratio-f unction  at 
a^j  is  directly  proportional  to  tan 
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The  assumed  properties  of  Mq  imply  the  following: 

(i)(R.(I)  7/  is  a non-empty  proper  subset  of  B 

which  is  in  the  immediate  vicinity  of  the  bottom  and  the 
right-hand  boundary  of  R. 

(ii)  There  is  a proper  subset  SCB  consisting  of 

I 

points  of  discontinuity  of  the  approximate  form  as  shown 
in  Fig.  7.1,  with  a "center  of  singularity"  at  the  middle. 

(iii)  The  remainder  of  B is  a region  where  the 
first  property  of  Mq  is  satisfied. 

For  illustrative  purposes,  reduce  the  two-dimensional 
problem  to  a one-dimensional  one.  Thus  the  set  A is 
condensed  into  a set  A'  of  equally-spaced  points  along  the 
vertical  center  line  R (designated  by  R' ) , and  the  set  B 
is  similarly  condensed  into  a set  B'  along  the  median  of  T 
drawn  from  the  left-hand  side  vertex  of  T (designated  by  T'). 
The  mapping  Mq  implies  that  the  upper  p.oint  in  A'  goes  to 
the  far  left-hand  side  of  T' , the  next  point  in  A'  goes  to 
the  next  location  on  T'  and  so  on,  where  the  spacing  between 
the  points  along  the  median  is  determined  by  the  equality  of 
the  areas  of  the  trapezoids  associated  with  the  points  of  B'. 
Let  a'  be  any  point  in  A'  and  let  k(a')  define  the  (one- 


U Range  of  I. 
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Figure  7.1  Geometrical  representation  of  a two-dimensional  grain-pile 
transformation  of  the  infinite  bin. 
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dimensional)  ratio-function  at  a'.  Joining  the  pairs 
(a*,  Mo(a'))  by  a family  of  mutually  non-interesting 
catenaries,  it  is  observed  that  tan  a'  is  monotonically 
decreasing  versus  the  depth,  until  the  point  of  inter- 
section (s)  between  R'  and  T'  is  reached.  At  this  point, 
tan  a'  is  undefined.  The  decrease  of  tan  a'  is  emphasized 
as  the  point  s is  approached  from  above,  starting  at  a 
height  of  approximately  2 h above  the  bottom  (where  h is 

b b 

the  height  of  the  center  of  singularity  s,  above  the  bottom) . 

Since  k(a')  oc  tan  a',  k(y)  has  the  same  qualitative 
behavior  as  described  for  tan  a' . This  result  establishes 
the  decreasing  nature  of  the  dynamic  k(y)  with  depth. 

(2)  The  location  of  the  center  of  singularity. 

By  geometrical  considerations 


v 


2 Hp  tan  0 


2 2 
D H tan 


e - 


tan  0 


D 

8 


tan 


0 


(a)  for  the 

infinite  bin 


(7.3) 

(b)  for  a circular 
cylindrical  bin 


No  comparison  between  the  infinite  bin  and  practical 
cases  can  be  made;  however,  a comparison  with  the  circular 
cylindrical  bin  is  in  order. 


139 


Consider  a typical  circular  cylindrical  grain  bin: 


e = 30 


o 


tan  0 = 0.577. 


Insert  these  values  in  Eq.  (7.3a)  to  obtain 


h 


0.127  H. 


(7.4) 


That  is,  the  lateral  pressure  decreases  rapidly  toward 
an  elevation  which  is  about  12%  above  the  bottom.  The 
maximum  pressure  occurs  above  this  elevation  along  a 
region  whose  length  is  1%  to  2 h , or  0.18H  to  0.24H  which 
is  altogether  0.30  H to  0.36  H,  or  one-third  of  the  total 
height  above  the  bottom. 

7.2  TOWARD  COMPUTER  IMPLEMENTATION  OF  GPT 


became  increasingly  more  difficult  to  preserve  the  applied 
aspects  of  the  theory.  Therefore,  a search  was  started  for 
ways  of  preserving  both  the  applied  aspects  and  the  theo- 
retical concepts  of  the  GPT,  while  by-passing  the  abstract 


of  numerical  analysis  and  computer  algorithms.  It  was  found 


computer  techniques,  developed  recently  at  the  Los  Alamos 


Formulating  this  theory  in  topological  terms  it 


notation.  The  direction  to  look  into  was  chosen  in  the  field 


that  certain  numerical  techniques  coupled  with  some  unique 
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Scientific  Laboratories,  can  be  useful  toward  this  end.  The 
method  can  be  described  roughly  as  a hybrid  between  digital 
simulation  techniques  and  the  more  traditional  finite-difference 
techniques.  It  can  carry  out  calculations  of  non-steady  phe- 
nomena and  output  the  results  as  time-dependent  realistic 
motion  pictures. 

In  relation  to  this  problem,  the  actual  mechanism  of 
GPT  can  be  computed  and  the  transformation  displayed  as  a motion 
picture.  By-products  of  such  simulations  would  be  various 
models  of  dynamic  ratio-functions. 

A complete  library  (51  programs;  more  than  6000  IBM 
cards)  of  a Stromberg-Carlson  4020  Micrifilm  Recorder  was 
acquired,  edited,  compiled,  debugged,  tested,  and  made  com- 
patible with  an  IBM  360/50  HASP  system.  The  combined  system 
has  the  capability  to  simulate  GPT  and  produce  the  results 
as  16  mm  motion  pictures. 


141 


CHAPTER  8 


RECOMMENDATIONS  FOR  FURTHER  STUDY 

8.1  INTRODUCTION:  COMPUTER  SIMULATION  VERSUS  EXPERIMENTATION 
Computer  simulation  has  come  into  increasingly  wide- 
spread use  to  study  the  behavior  of  complex  systems  whose  state 
changes  over  time.  Alternatives  to  the  use  of  simulation  are 
mathematical  analysis,  experimentation  with  either  the  actual 
system  or  a prototype  of  the  actual  system,  or  reliance  upon 
experience  and  intuition.  All,  including  simulation,  have 
limitations.  Mathematical  analysis  of  complex  systems  such 
as  the  grain-bin-system  is  very  often  impossible;  experimenta- 
tion with  actual  or  pilot  systems  is  costly  and  time  consuming, 
and  the  relevant  parameters  are  not  always  subject  to  control. 
Intuition  and  experience  are  often  the  only  alternatives 
available,  but  can  be  very  inadequate. 

Simulation  problems, are  characterized  by  being  mathe- 
matically intractable  and  having  resisted  solution  by  analytic 
methods.  The  problems  usually  involve  many  variables,  many 
parameters,  functions  which  are  not  well-behaved  mathematically, 
and  random  variables.  The  simulation  is  often  a technique  of 
last  resort.  Yet,  much  effort  is  now  devoted  to  computer 
simulation  because  it  is  a technique  that  gives  answers  in 
spite  of  its  difficulties,  cost,  and  time. 

During  this  study  various  computer  simulation  techniques 
were  developed  and  exploited  with  great  success.  However, 
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at  various  points  (for  example,  bottom  of  page  44)  it  was 
pointed  out  that  ultimately  empirical  determination  of  certain 
parameters  must  be  made  and  then  fed  into  the  basic  model 
in  order  to  establish  practical  solutions.  In  fact,  it  is 
desirable  that  both  simulation  "experiments”  and  real  experi- 
ments be  carried  out,  the  results  compared,  the  techniques 
adjusted,  etc.,  until  satisfactory  solutions  are  obtained. 

This  procedure  is  summarized  in  the  following  diagram. 
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As  can  be  seen  in  the  diagram,  a dual  computer  and 


experimentation  approach  is  proposed  for  the  study  of  grain 
pressures.  On  the  one  hand  various  computer  simulation 
experiments  are  carried  out,  involving  choices  of  geometries, 
characteristic  functions  and  other  parameters.  On  the  other 
hand,  actual  experimentation  is  carried  out,  also  involving 
various  conditions.  The  results  of  the  computer  analyses 
and  the  actual  experiments  are  compared,  and  the  differences 
between  predicted  and  observed  results  are  noted.  These,  in 
turn,  indicate  changes  to  be  made  in  the  computer  experiments 
as  well  as  the  real  experiments,  and  so  forth. 

The  desired  goal  is,  of  course,  a self-consistent 
theory  encompassing  all  known  factors  affecting  grain  pressures. 
The  computer  simulation  and  the  engineering  experimentation 
results  are  integrated,  synthesized,  and  tested  in  an  overall 
computer  simulation  of  grain  pressures.  The  simulation  has 
the  advantage  of  demonstrating  the  total  relation  of  many 
and  varied  factors  that  contribute  to  the  overall  grain  pressure 
mechanism;  it  enables  the  quantitative  study  of  the  relation- 
ship among  many  parts  and  factors,  and  the  evaluation  of  the 
contribution  of  any  single  factor  to  the  total  picture;  in 
addition,  the  simulation  enables  the  study  of  the  dynamic  or 
changing  situation  which  is  intrinsic  to  grain  pressure  systems. 
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8.2 


RECOMMENDATIONS  FOR  EXPERIMENTAL  INVESTIGATIONS 


8.2.1  Grain  Mechanics 

a.  General . Conduct  laboratory  investigations  on 
the  geometric,  mechanical  and  rheological  properties  of 
grain,  by  methods  analogous  to  those  of  soil  mechanics, 
rheology  and  related  fields.  Also,  investigate  the  special 
properties  of  grain  due  to  organic  and  biological  factors, 
as  affecting  grain  pressures  in  storage  structures. 

b.  Specific.  The  following  is  a list  of  some  grain 
properties  that  were  found  to  be  of  some  significance  in 
studying  grain  pressure  effects: 

The  single  kernel:  Average  weight,  average  volume,  specific 
weight,  average  surface  area,  typical  geometric  form,  surface 
roughness,  hardness  (friction  resistance) , breakage  and 
mechanical  damage,  other  irregularities,  modulus  of  elasticity, 
compressive  strength,  shear  strength,  impact  strength. 

Bulk  of  grain:  Unit  density,  void  ratio  and  porosity,  angle 
of  natural  slope,  internal  friction,  coefficients  of  friction 
on  various  materials,  moisture  content  (under  various  tem- 
peratures and  humidities) , modulus  of  elasticity,  compressive 
strength,  static  and  dynamic  shear  strength,  modulus  of 
resilience,  modulus  of  toughness,  stress  relaxation,  flow- 
factor,  susceptibility  to  vibration,  consolidation,  specific 
heat,  diaelectric  constant,  thermal  conductivity,  coefficient 
of  volume  expansion,  organic  and  biological  variations  in 
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volume,  other  biological  properties  (such  as:  absorption  of 
moisture,  generation  of  heat,  carbon  dioxide,  etc.),  per- 
centage of  foreign  matter  and  dust. 

8.2.2  Full-Scale  Test  Structures 

It  is  recommended  to  design  and  construct  test  bins 
to  carry  out  series  of  controlled  experiments.  A typical 
test  bin  should  be  of  circular  cross-section,  built  of 
reinforced  concrete  or  steel,  and  have  sufficiently  large 
dimensions,  i.e.,  at  least  15  ft.  diameter  and  75  ft.  height. 
Provisions  for  changing  the  hopper  should  be  incorporated 
into  the  design.  At  a minimum,  the  test  bin  should  be 
instrumented  for  measuring  the  following  parameters: 

(1)  lateral  pressure 

(2)  vertical  load  on  the  walls 

(3)  vertical  and  lateral  pressure  within  the  grain 

(4)  bottom  pressure 

(5)  grain  density 

The  measuring  devices  should  be  of  adequate  range  and 
sensitivity  and  suitable  to  record  time-dependent  variations 
in  pressures.  Their  number  and  distribution  should  be  adequate 
to  cover  the  entire  system  and  to  provide  good  resolution, 
capable  of  matching  the  resolution  of  the  computer  simulation. 

The  results  of  the  measurements  should  be  used  to 
determine  empirically  the  characteristic  functions  under 
various  conditions,  static  and  dynamic.  In  particular, 
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a series  of  experiments  should  be  conducted  to  satisfy  the 
following  specific  objectives: 

a.  Charging  - To  determine  the  relation  of: 

(1)  The  point  of  charging  - center  or  side 

(2)  Method  of  charging  - continuous  or  intermittent 

(3)  Rate  of  charging  - various  levels,  constant  or 
varying 

(4)  "Rain  charging"  - various  methods 

(5)  Inclined  flow  charging  through  the  wall 

(6)  Charging  by  special  mechanical  devices  to  reduce 
kinetic  energy.  For  example:  perforated  pipe 
according  to  the  Reimberts'  method, 

to  the  vertical  loads  and  lateral  and  bottom  pressures, 

in  the  static  and  dynamic  states. 

b.  Discharging  - To  determine  the  relation  of: 

(1)  The  point  of  discharging  - central,  to  the  side, 
at  the  wall,  many  gates  at  various  distributions 

(2)  Various  hopper  slopes 

(3)  Method  of  discharging  - continuous  or  intermittent 

(4)  Rate  of  discharging  - various  levels,  constant  or 
varying 

(5)  Perforated  discharging  pipe  - central  or  side 

(6)  Horizontal  rings  along  the  perimeter  - various 
dimensions  and  spacings 

(7)  Various  conditions  of  gate  shut-off  time. 
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to  the  vertical  loads  and  lateral  and  bottom 
pressures,  in  the  static  and  dynamic  states. 

c.  Flow  Pattern  - To  determine  the  relation  of  the 
flow  pattern  during  charging  and  discharging  to  the 
above  mentioned  loads  and  pressures. 

d.  Simultaneous  Charging  and  Discharging  - To  determine 
the  relation  of  simultaneous  charging  and  discharging 
to  the  above  mentioned  loads  and  pressures. 

e.  Grains  - To  determine  the  relation  of: 

(1)  Various  types  of  grain 

(2)  Various  densities 

(3)  Coefficient  of  friction  of  grain  on  grain 

(4)  Moisture  content 

(5)  Cracked  grain  and  foreign  material, 

to  the  above  mentioned  loads  and  pressures. 

f.  Aeration  - To  determine  the  relation  of: 

(1)  Grain  temperature  (changes  caused  during  cooling 
and  warming) 

(2)  Small  changes  in  grain  moisture  content  (both  for 
increasing  and  decreasing  the  moisture  content) , 

to  the  vertical  loads  and  lateral  and  bottom  pressures 
in  the  static  state. 

g.  Vibration  - To  determine  the  effect  of  forced  vibration 
producing  settling  on  the  vertical  loads  and  lateral 
and  bottom  pressures,  in  the  static  and  dynamic  states. 
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h.  Wall  Surface  - To  determine  the  effect  of  various 
wall  surfaces: 

(1)  Smooth  steel 

(2)  Smooth  and  rough  coating  on  steel 

(3)  Concrete, 

on  the  vertical  loads  and  lateral  and  bottom  pressures, 
in  the  static  and  dynamic  states. 

i.  Period  of  Storage  - To  determine  the  effect  of  long 
storage  periods  (a  year  or  longer) , on  the  vertical 
loads  and  lateral  and  bottom  pressures  in  the  static 
state . 

j.  Solar  Radiation  - To  determine  the  effect  of  solar 
heating  on  the  vertical  loads  and  lateral  and  bottom 
pressures  in  the  static  state. 
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APPENDIX  B 


PROGRAM  LISTINGS  FOR  THREE-DIMENSIONAL  PLOTS 

Under  this  project,  a computer-based  automatic  system 
to  plot  any  two-variable  function  as  a three-dimensional 
surface  was  implemented.  The  system  is  an  adaptation  from 
a NASA  Technical  Memorandum  (TM  x-1598,  June  1968).  It. 
utilizes  an  IBM  360/50  computer  and  a Calcomp  digital  incre- 
mental plotter.  The  software  is  written  in  Fortran  IV. 

In  section  5.2  various  analytic  models  are  displayed 
graphically  as  three-dimensional  surfaces.  Figures  5.1  through 
5.12  are  examples  of  the  capabilities  of  the  system.  The  same 
system  can  be  used  to  describe  lateral  pressure  diagrams  as 
three-dimensional  surfaces  over  the  walls,  two-dimensional 
characteristic  functions,  etc. 

The  program  listings  follow.  In  addition,  the  following 
standard  Calcomp  subroutines  are  assumed  to  be  available: 

PLT770 , NUMBER,  SYMBOL. 
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THREE-DIMENSIONAL  PLOTS  PACKAGE 


EXTERNAL  FB  ... 

LOGICAL  CUBE  » DR A WMF 
DIMENSION  S! 6000) , IBUFF(6000) 

CALL  PLOTS ( I BUFF  ,6000, 1 ) 

CALL  PLOT ! 0*  0,— 30®0,3) 

CALL  PLOT (0. 5,-29® 5,-3 ) 

CUBE  = *TRUE*  _ 

CALL  PL0T3D(0o  , 120o ,0. ,400o  , S , 51 , 31, 7, 51 , FB, CUBE ) 
DRAWME  = oFALSEo 

CALL  ROTATE! 0*0 ,0. 0*45.0 , DRAWME ) 

DRAWME  = •TRUE, 

DO  2 K=  1 , 3 

CALL  ROTATE!0®0, 20*0, 0.0, DRAWME) 

NBL=NBL+1 
2 CONTINUE 

CALL  WHERE! XTEM,YTEMtDUM) 

CALL  PLOT (XTEMt YTEM,999) 

STOP 

END 


FUNCTION  FB ! X , Y ) 

DATA  GAMAO/45® /, DELTA /. 06/, FKO/* 333/, FMU/. 425/, R/5./, FLAMDA/. 425/ 

GAMMA=G AMAO  + DELTA*X 

FK  = 2**FK0  / !la+  ! FMU*FKO*X/R ) ) 

FB  = GAMMA  - FLAMDA  * FK  * Y / R 

RETURN 

END 


BLOCK  DATA 

COMMON  /XCPIDX/  NBL , RUNNO ,CP I D ( B ) 

DATA  NBL /I / , RUNN0/4HNQ®  9/ , CPID/4HISTR,7*4H 
END 
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_ SUBROUTINE  PLOT3D ( X M I N , XMAX , YMI N , YMAX , S , NXPTS , NYPLS , 
* NXPLS,  NYPTS  » F » CUBE ) 


PLOT3D 

PLOT3D 


COMMON  / LFACT/ 
LOGICAL  CUBE 
DIMENSION  S( 1 ) 
DATA  IK/O/ 
IFJIK.NEoOJ  GO 
CALL  WRITE 
I K=  1 
FNX1 
FNY1 
FNX2 
FNY2 
LI  = 

L2  = 

N1  = L1 

N2=N1 +L 1 _ 
N3=N2+L1 
N4=N3+l 2 
N5=N4+ 12 
DX=XMAX“XMIN 
DY= YMAX-YMI N 


FACT 


TO  1 


= NXPTS-1 
= NYPLS-1 
= NXPLS”! 

= NYPTS-1 

NXPTS*NYPLS 

NXPLS*NYPTS 


+ 1 


DO  2 1 = 1 » NXPTS 
002  J=l, NYPLS 
II  = NXPTS*«J-1  ) 

11  = 1 I 

I 2 = Ni  + 1 I 

I 3=N2+  1 1 

_SU  1 )=XMIN«-DX*FLOAT(  1-1 ) /FNX1 
S(I2)=YMIN+DY*FL0AT( J-lJ/FNYl 
2 S ( I 3 ) = F( S ( 1 1 ) * S ( 1 2 ) I 


DO  3 1=1, NXPLS 
DO  3 J=l, NYPTS 
II  = NXPLS* ( J-l ) 


+ I 


PL0T3D 
PL0T3D 
PL0T3D 
PL0T3D 
PL0T3D 
PL0T3D 
PL0T3D 
PL0T3D 1 0 
PL0T3D1 1 
PL0T3D12 
PL0T3D1 3 
PL0T3D14 
PL0T3D1 5 
PL0T3D16 
PL0T3D17 
PL0T3D1 8 
PL0T3D1 9 
PLOT3D2  0 
PL0T3D2 1 
PLOT3D22 
PLOT3D23 
PL0T3D24 
PL0T3D2  5 
PL0T3D26 
PLOT3D27 
PL0T3D2  8 
PL0T3D29 
PL0T3D  30 
PLOT3D3 1 
PLOT3D32 
PL0T3D33 
PL0T3D34 
PL0T3D35 
PL0T3D  3 6 
PLOT3D3  7 


I1=N3«-II 

PL0T3D3  8 

I 2=N4+ I T 

PL0T3D39 

1 3=N5+ 1 1 

PL0T3D40 

SII  1)  = XMIN4-DX*FL0AT(I-1)  /FNX2 

PL0T3D41 

S < I 2 ) = YMI N+D Y*FL0AT ( J”1 ) / FNY2 

PL0T3D42 

3 S ( I 3 ) =F( S ( 1 1 ) , S ( I 2 ) ) 

PLOT  3D43 

PL0T3D44 

PLOT  3D45 

N1=N1+1 

PL0T3D46 

N2=N2+ 1 

PL0T3D47 
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N3=N3+1 

N4=N4+1 

N5=N5+1 

CALL  SCALE ( XMIN, XMAX , YMI N , YMAX , S( 1 ) f S(N1 ) , S ( N2 ) , NX PTS , NYPL S , 
* S(N3) , S ( N4 ) * S ( N5 ) , NXPLS ,NYPTS , CUBE ) 

C 

C 

CALL  AXIS(0, . FALSE. ) 

C 

C 

A SUM=0»  0 .... 

9 B$UM=0» 

. C SUM=0®  . 

RETURN 

C 

s C 

ENTRY  ROTATE! ALPHA, BETA, GAMMA, DRAWME) 

LOGICAL  DRAWME 

IF  (FACT  .NE.  lo ) CALL  FACTOR  (FACT) 

CALL  PLOT (Oo  0 » 0,0,  -3) 

CALL  PLOT (0,0,  5,0,  -3) 

IF  (FACT  .NE.  I, ) CALL  FACTOR  (1.) 

CALL  TRNMAT  ( ALPHA , BET A , GAMMA ) 

CALL  PHI  (S(1),S(N1)  ,S(N2) ,NXPTS,NYPLS) 

CALL  PHI  ( S ( N3 ) , S ( N4 ) ,S(N5) , NX PL S, NY PTS I 

CALL  AXIS! 1, DRAWME) 

ASUM=ASUM«-ALPHA 

BSUM=BSUM+BETA 

CSUM=C  SUM+GAMMA 

IF( 9 NOT # DRAWME ) RETURN 

CALL  DRAW$(S(1),S(N1) , S ( N2 ) , NXPTS , NYPL S ) 

CALL  DRAW(S(N3) ,S(N4) ,S(N5) , NXPLS , NY PTS ) 


CALL  WRITES( ASUM,BSUM,CSUM) 

IF  (FACT  oNE,  lo ) CALL  FACTOR  (FACT) 
CALL  PLOT (10o0»-5o0,-3) 

IF  (FACT  .NE.  I,)  CALL  FACTOR  (1.) 

RETURN 

END 


PL0T3D48 
PL0T3D49 
PLOT3D  50 
PLOT3D5 1 
PL0T305  2 
PLOT3D53 
PLOT3D54 
PLOT3D55 
PL0T3D56 
PLOT3D57 
PL0T3D5  8 
PLOT3D5  9 
PLOT3D60 
PL0T3D61 
PLOT3D62 
PLOT3D63 
PL0T3D64 
PL0T3D6  5 

PL0T3D66 

PL0T3D67 

PL0T3D68 
PL0T3D6  9 
PL0T3D  70 
PLOT3D71 
PLOT3D72 
PL0T3D73 
PL0T3D74 
PL0T3D75 
PLOT3D76 
PL0T3D77 
PL0T3D78 
PL0T3D7  9 
PL0T3D  80 

PLOT3D81 

PLOT3D82 

PLOT3D83 


BLOCK  DATA  1 

COMMON/LABEL/  LAB(18)  2 

COMMON  /LISTR/  ISTR 
COMMON  /LFACT / FACT 

DATA  ISTR  /4HISTR/  __  _ 

DATA  FACT/1. / 

DATA  LAB  /4HISTR,  17*4H 

END  * 
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SUBROUTINE  DRAW  ( X , Y , Z , NX , NY ) 

COMMON  /LFACT/  FACT 

IF  (FACT  • NE o 1®)  CALL  FACTOR 

(FACT) 

DRAW 

1 

DIMENSION  X ( NX , NY ) , Y ( NX , NY ) , 

Z ( NX , NY ) 

DRAW 

2 

00  1 1=1, NX 

DRAW 

3 

CALL  PLOT  ( Y ( I ,1) ,Z( 1,1), 3) 

DRAW 

4 

DO  1 J=2 , NY 

DRAW 

5 

CALL  PLOT  (Y( I , J)  ,Z( I , J) ,2) 

DRAW 

6 

CONTINUE 

DRAW 

7 

IF  (FACT  ® NE®  1®  ) CALL  FACTOR 
RETURN 

(1®  ) 

DRAW 

8 

ENTRY  DRAWS  ( X , Y , Z , NX , NY ) 

DRAW 

9 

IF  (FACT  ® NE  ® 1®  ) CALL  FACTOR 
DO  2 J = 1 » NY 

(FACT) 

DRAW 

10 

CALL  PLOT  (Y(ltJ)fZ(lfJ) ,3) 

DRAW 

11 

DO  2 1=2, NX 

DRAW 

12 

CALL  PLOT  (Y( I , j) ,Z( I , J) ,2) 

DRAW 

13 

CONTINUE 

DRAW 

14 

IF  (FACT  oNE,  1®)  CALL  FACTOR 
RETURN 

( 1®  ) 

DRAW 

15 

END 

DRAW 

16 

SUBROUTINE  WRITE 

WRITE 

1 

COMMON  /LFACTZ  FACT 
COMMON  /LISTR/  ISTR 
COMMON  /MAXES/  XMIN, XMAX,YMIN 

,YMAX,ZMIN, ZMAX 

WRITE 

2 

COMMON  /XCPIDX/  NBL , RUNNO, CPID( 8 ) 

WRITE 

3 

COMMON  /LABEL/  LAB(18) 

WRITE 

4 

DIMENSION  BLOCK (2) , AL PHA ( 2 ) , BE T A ( 2 ) , GAMM A ( 2 > 

,HXMIN( 2) ,HXMAX( 2) , WRITE 

5 

HYMIN( 2) ,HYMAX( 2) , H ZM  I N ( 2 ) , HZM AX ( 2 ) 

WRITE 

6 

DATA  BLOCK, ALPHA, BETA, GAMMA/4HBL0C , 4HK  ,4HALPH,4HA=  , 4HBETA , 4H=WR I TE 

7 

, 4HG AMM  ,4HA=  / 

WRITE 

A7 

DATA  HXMIN,HXMAX ,HYMI N, HYMAX, 

HZMIN»H  ZM AX/4HXMIN ,4H=  ,4HXMAX,4H=  WRITE 

B7 

, 4HYM IN, 4H=  ,4HYMAX,4H= 

, 4H  Z MI N,4H= 

,4HZMAX,4H=  WRITE 

C7 

IF  (CPID(l). EQo ISTR)  RETURN 
IF  (FACT  oNE®  1®)  CALL  FACTOR  (FACT) 
CALL  PLOT  ( 0®  0 ,0a  0 ,-3  ) 

CALL  PLOT  (OoO,5oO,-3) 

C P I D ( 8 ) = RUNNO 

CALL  SYMBOL  ( Oo 0 , -4. 0 , ® 25 , C P I D , 90 ® 0, 32 ) 
CALL  SYMBOL  (0®0,-5o0,0o07,BL0CK,0«0,6) 
NBL1=NBL-1 
FNB  L=NBL 1 

CALL  NUMBER(0®63,-5e0,0.07,FNBL,0®0,-l ) 

WRITE  (6,2)  NRLI 

CALL  PLOT  ( 9®  0 » -5o  0 » -3  ) 

IF  (FACT  ® NE®  1®)  CALL  FACTOR  (1®) 
RETURN 


WRITE 

WRITE 

WRITE 

WRITE 

WRITE 

WRITE 

WRITE 

WRITE 

WRITE 

WRITE 


8 

9 

10 

11 

12 

13 

14 

15 

16 
17 


WRITE  18 
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WRITE  19 


ENTRY  WRITES  (A1,B1,C1) 

IF  (FACT  .NE.  lo ) CALL  FACTOR  (FACT) 

A=AMOD (Al, 360.0)  ...  ...  WRITE  20 

B=AM0D( B1 ,360.0 > WRITE  21 

C=AM0D(C1 ,360.0)  WRITE  22 


CALL 

SYMBOL 

(0.0, -5. 0,0. 07, BLOCK ,0.0,6) 

WRITE 

23 

FN8L  = 

:NBL 

WRITE 

24 

CALL 

NUMBER (0,6 3,-5. 0,0. 07,FNBL ,0. 0,-1 ) 

WRITE 

25 

WRITE 

(6,2) 

NBL 

WRITE 

26 

IF  (LAB( 1 ). EQo ISTR ) GO  TO  1 

WRITE 

27 

CALL 

SYMBOL 

(-3. 0,-4. 0,0.1 , LAB, 0.0, 72) 

WRITE 

28 

CALL 

SYMBOL 

( -3. 0, -4. 2 5, 0. 1, ALPHA, 0.0, 6) 

WRITE 

29 

CALL 

NUMBER 

(— 2.4,— 4.25,0.1,A,0o0,5) 

WRITE 

30 

CALL 

SYMBOL 

(-3. 0,-4. 5,0. 1, BET A, 0.0, 6) 

WRITE 

31 

CALL 

NUMBER 

(-2, 4, -4.5  ,0. 1 , B,0. 0, 5) 

WRITE 

32 

CALL 

SYMBOL 

(-3.0, -4075, 0.1, GAMMA, 0. 0,6) 

WRITE 

33 

...  CALL 

NUMBER 

(-2. 4, -4.75,0.1,0,0.0, 5) 

WRITE 

34 

CALL 

SYMBOL 

(-1,0,— 4. 25,0.1,HXMIN,0.0,5) 

WRITE 

35 

CALL 

NUMRER 

(-0.4,-4.25,0.1,XMIN,0.0,5) 

WRITE 

36 

CALL 

SYMBOL 

(-1.0, -4.5  ,0.1,HYMIN,0.0,5) 

WRITE 

37 

CALL 

NUMBER 

(-0.4, -4.5  ,0. 1, YMIN,0.0, 5) 

WRITE 

38 

CALL 

SYMBOL 

( — loO,— 4.75, 0.1»HZMIN, 0.0,5) 

WRITE 

39 

CALL 

NUMBFR 

(-0, 4,-4. 75,0. 1,ZMIN,0.0,5) 

WRITE 

40 

CALL 

SYMBOL 

( 1.0,-4.25,0.1,HXMAX,0.0,5) 

WRITE 

41 

CALL 

NUMBER 

( 1.6, -4.25,0. 1 , XMAX,0.0,5) 

WRITE 

42 

CALL 

SYMBOL 

( 1.0, -4. 5 , 0. 1 , HYM AX ,0.0,5) 

WRITE 

43 

CALL 

NUMBER 

( 1.6, -4.5  ,0. 1,YMAX ,0.0,5) 

WRITE 

44 

CALL 

SYMBOL 

( 1.0,— 4.75,0. 1 ,HZMAX,0. 0, 5 ) 

WRITE 

45 

CALL 

NUMBER 

( 1.6, -4*75, 0.1, ZMAX, 0.0,5) 

WRITE 

46 

IF  (FACT  .NE 

. 1.  ) CALL  FACTOR  ( 1.  ) 

RETURN 

WRITE 

47 

WRITE 

48 

2 . FORMAT  ( 1 20  X 

,5HBL0CK, 15 ) 

WRITE 

49 

END 

WRITE 

50 
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SUBROUTINE  SCALE  ( XM I N , XMAX , YM I N , YMA X , X 1 , Y1 , Z 1 , NX 1 f N Y1 , X 2 » Y2 , Z2 , NXSC ALE  1 


*2»NY2,CURE) 

SCALE 

2 

DIMENSION  Xl(NXIfNYl)t  Y1(NX1,NY1),  Z1(NX1,NY1), 

X2 ( NX2 » NY2 ) , 

Y2 (NSCALE 

3 

*X2,NY2>,  Z2 ( NX2 1 NY2 ) 

SCALE 

4 

COMMON  /MAXES/  XM I N1 , X MAXI , Y MI N I ♦ Y MAXI t ZM I N 9 ZMAX 

SCALE 

5 

LOGICAL  CUBE 

SCALE 

6 

REAL  MAXXtMAXY, MAXZ 

SCALE 

7 

XMIN1=XMI N 

SCALE 

8 

XMAX1=XMAX 

SCALE 

9 

YMI N1=YMIN 

SCALE 

10 

YMAX1=YMAX 

SCALE 

11 

MAXX=(XMAX-XMIN)/3o5 

SCALE 

12 

MAXY=( YMAX-YMIN ) / 3*  5 

SCALE 

13 

ZMAX=Z1 (1,11 

SCALE 

14 

ZMI N=Z 1 ( 1,1) 

SCALE 

15 

DO  1 1=1, NX1 

SCALE 

16 

DO  1 J=! , NY I 

SCALE 

17 

ZMA  X = AMAX1 ( Zl(  I » J)  ,ZMAX) 

SCALE 

18 

1 

ZMIN=AMIN1 (Z1 ( I , J) , ZMIN) 

SCALE 

19 

DO  2 1=1, NX2 

SCALE 

20 

DO  2 J=1,NY2 

SCALE 

21 

ZMAX=AMAX1 (Z2(T,J),ZMAX) 

SCALE 

22 

2 

ZMI N= AMI N1 ( Z2 ( I, J), ZMIN) 

SCALE 

23 

MAXZ=( ZMAX-ZMIN) /3, 5 

SCALE 

24 

IF  (CUBE)  GO  TO  3 

scale 

25 

MAX X= AM A XI (MAXX»MAXYtMAXZ) 

SCALE 

26 

MAXY=MAXX 

SCALE 

27 

MAXZ=MAXX 

SCALE 

28 

3 

DO  4 1=1, NX1 

SCALE 

29 

DO  4 J=1 , NY1 

SCALE 

30 

X1(I,J)-(X1( I »J )-( XMAX+XMIN) /2o )/MAXX 

SCALE 

31 

Y1 ( I , J)= ( Y1 ( I T J )-( YMAX+YMIN) /2, ) /MAXY 

SCALE 

32 

4 

Z1(1,J)=(Z1< I tJ)~( ZMAX+ZMIN) /2« ) /MAXZ 

SCALE 

33 

DO  5 1=1, NX2 

SCALE 

34 

DO  5 J= 1 , NY2 

SCALE 

35 

X2(I ,J)=(X2(I, J)-( XMAX+XMIN) /2, ) /MAXX 

SCALE 

36 

Y2(I , J)=(Y2( I , J)-( YMAX+YMIN>/2, )/MAXY 

SCALE 

37 

5 

Z2(I»J)=(Z2(I , J)-( ZMAX+ZMIN)/2, ) /MAXZ 

SCALE 

38 

RETURN 

SCALE 

39 

END 

SCALE 

40 

172 


SUBROUTINE  AX  I S ( I f DRAW ME ) 

AXIS 

1 

COMMON  /L F ACT/  FACT 

LOGICAL  DRAWME 

AX I S 

2 

DIMENSION  X(4,l)  ,Y{4,l)*Z(4f 1) 

AXIS 

3 

IF  ( I.NE.O)  GO  TO  1 

AXIS 

4 

X(ltl]*0, 

AXIS 

5 

X(2» 1 )=loO 

AXIS 

6 

X ( 3 * 1 ) = 0. 

AXIS 

7 

X (4» 1 1 =0® 

AXI  S 

8 

Y ( 1 t 1 ) =0. 

AXIS 

9 

Y ( 2 » 1 ) =0o 

AXIS 

10 

> 

Y(3»l)=l»0 

AXIS 

11 

Y ( 4» 1 ) =0o 

AXIS 

12 

Z ( 1 » 1 ) =0o 

AXIS 

13 

Z ( 2 » 1 > =0o 

AXIS 

14 

Z(3,,l  )=0. 

AXIS 

15 

Zr*f » 1 ) -1.0 

AXIS 

16 

RETURN 

AXIS 

17 

l 

CALL  PHI (X,YfZ,4,l ) 

AXI  S 

18 

IF  (DRAWME)  GO  TO  2 

AXIS 

19 

RETURN 

AXIS 

20 

2 

YY=Y(1  ,1  ) 4-4.0 

AXIS 

21 

ZZ=Z  ( 1 1 1 ) -t-4*  0 

AXIS 

22 

IF  (FACT  .NE.  1.)  CALL  FACTOR  (FACT) 

CALL  PLOT ( YY  * Z Z t 3 ) 

AXIS 

23 

CALL  PL0T(Y(2,1 )+4. , Z(2,l)+4»,  2) 

AXIS 

24 

CALL  SYMB0L(Y(2 ,1 )+3.  97143,  Z(2,l)*3.95, 

• It 

1HX, 

o 

• 

o 

1 ) 

AXIS 

25 

CALL  PL0T(YYtZZ,3) 

AXIS 

26 

CALL  P L0T(Y(3,1 )+4. , Z(3,l)+4o,  2) 

AXIS 

27 

CALL  SYMROL  (Y  ( 3 » 1 ) 3<»  97 1 43 » Z(3,l)+3.95f 

• It 

1HY, 

o 

. 

o 

1) 

AXI  S 

28 

CALL  PLOT ( YY  t Z Z * 3 ) 

AXIS 

29 

CALL  PL  OT ( Y ( 4 » 1 ) +4,  , Z(4,l)+4.,  2) 

AXIS 

30 

CALL  SYMBOL  (Y(4»  1 )*3.97143f  Z(4fl)+3.95, 

• If 

1HZ  f 

o 

. 

o 

«• 

1 ) 

AXIS 

31 

IF  (FACT  .NE.  U ) CALL  FACTOR  (1.) 

RETURN 

AXIS 

32 

END 

AXIS 

33 
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SUBROUTINE  TRNMAT  ( AL PHA , B ET A , GAMMA ) 

COMMON  /MATRIX/  TMAT(3,3) 

A=ALPHA/57„2957795 

B=  BET A/57» 2957795 

C=GAMMA/570 2 957795 

S INA=S IN( A ) 

S I NB  = S I N ( B ) 

SI NC=S I N ( C ) 

COSA=COS( A) 

COSB=C  OS ( B ) 

COSC=COS(C) 

TMAT(1 ,l)=COSC*COSB 

TMAT(l,2)=COSC*SINB*SINA-SINC*COSA 
TMATC1 ,3) =COSC*S INB*COSA+S INC*S INA 
TMAT<2, 1 )=SINC*COSB 

TMAT(2,2)=SINC*SINB*SINA+C0SC*C0SA 

TMAT I 2,3) =SINC*SINB*COSA-COSC*S INA 

TMAT ( 3 , 1 ) =- S I NB 

TMAT(3,2  >=COSB*SINA 

TMAT(3,3)=COSB*COSA 

RETURN 

END 


TRNMAT01 
TRNMAT02 
TRNMAT  3 
TRNMAT  4 
TRNMAT  5 
TRNMAT  6 
TRNMAT  7 
TRNMAT  8 
TRNMAT  9 
TRNMATIO 
TRNMAT1 1 
TRNMAT12 
TRNMAT1 3 ; 

TRNMAT14  4 
TRNMAT1 5 
TRNMAT1 6 v 
TRN^AIJJ 
TRNMAf 18 
TRNMAT1 9 
TRNMAT20 
TRNMAT2 1 
TRNMAT22  ’ 


SUBROUTINE  PHI  ( X, Y, Z ,NX ,NY ) 

DIMENSION  X(NX,NY),  Y(NX,NY),  Z ( NX ,NY ) 


COMMON  /MATRIX/  TMAT (3, 3) 


DO  1 1=1, NX 
DO  1 J=1 , NY 


XP=TMAT ( 1 , 1 ) *X ( I »J)+TMAT(1,2)*Y(I»J 
YP=TMAT  ( 2 ♦ 1 ) *X  (I  , J)  +TMAT  (2» 2)*Y( I ♦ J 
ZP  = TMAT ( 3 , 1 ) *X ( I , J)+TMAT ( 3, 2)*Y(  I,  J 


X ( I , J)=XP 
Y ( I , J)=YP 
1 _Z(I , J)=ZP 


RETURN 

END 


174 


PHI 

3 

PHI 

4 

PHI 

5 

PHI 

6 

PHI 

7 

PHI 

8 

PHI 

9 

) +TMAT  ( 1 * 3 )*Z ( 

I 

, J I 

PHI 

10 

) +TMAT ( 2 , 3 )*Z ( 

I 

tJ) 

PHI 

11 

) +TM AT ( 3 , 3 )*Z ( 

I 

PHI 

12 

PHI 

13 

PHI 

14 

PHI 

15 

PHI 

16 

PHI 

17 

PHI 

18 

PHI 

19 

PHI 

20 

PHI 

21 

t 

t 


UNITED  STATES  DEPARTMENT  OF  AGRICULTURE 
AGRICULTURAL  RESEARCH  SERVICE 
H YATTSVILLE,  MARYLAND  20782 


OFFICIAL  BUSINESS 


POSTAGE  & FEES  PAID 
United  States  Department  of  Agriculture 


